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ABSTRACT 

The  problem  of  a  porous  graphite  fiber  plate  subject  to 
combustion  is  formulated.   The  transient  one-dimensional 
model  leads  to  two  heat  transfer  equations  on  the  graphite 
and  air,  and  a  mass  transfer  equation  on  the  oxygen.   In 
addition  to  a  combustion  heat  generation  term,  the  heat 
transfer  model  includes  the  mechanisms  of  conduction,  radia- 
tion between  fibers,  and  convection  due  to  induced  air  flow 
through  the  mat.   The  mass  transfer  model  includes  diffusion 
and  convection  mechanisms,  as  well  as  a  combustion  consumption 
term.   The  temperature  dependency  of  the  system  parameters 
is  taken  into  account.   The  resulting  nonlinear,  coupled 
partial  differential  equations  are  solved  by  a  Galerkin  formu- 
lation of  the  finite  element  method.   A  number  of  computer 
analysis  results  are  presented.   The  results  show  the  effects 
of  initial  conditions,  plate  thickness,  fiber  diameter  and 
air  flow  rate  on  ignition  and  extinction. 
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I.   INTRODUCTION 

The  combustion  of  a  composite  laminate  plate  consisting 
of  graphite  fibers  in  an  epoxy  matrix  is  discussed.   The 
combustion  process  for  a  graphite  laminate  takes  place  in 
two  stages,  epoxy  combustion  followed  by  combustion  of  the 
graphite.   At  the  end  of  the  first  stage  the  epoxy  has  burned 
away  exposing  a  porous  graphite  mat.   Spacing  between  the 
fibers  is  maintained  by  the  residue  of  the  combustion  pro- 
ducts.  This  present  work  is  concerned  with  the  last  stage 
of  the  combustion  process,  specifically,  the  thermal  response 
of  the  remaining  graphite  fiber  mat. 

The  selection  of  the  second  stage  for  the  analysis  was 
made  as  a  result  of  observations  during  composite  plate 
"burn"  experiments  at  the  Naval  Weapons  Center,  China  Lake  [1] 
The  objective  of  the  experiments  was  to  assess  damage  of  com- 
posite aircraft  structures  subjected  to  open  deck  fires  on 
aircraft  carriers.   For  the  experiments,  an  epoxy-graphite 
laminate  was  mounted  into  the  wall  of  a  wind  tunnel  parallel 
to  the  flow.   One  side  of  the  composite  plate  was  exposed  to 
the  wind  tunnel  flow  and  the  opposite  side  was  open  to  the 
environment.   This  arrangement  was  to  simulate  a  composite 
wing  section  or  fuselage  subjected  to  typical  wind  condi- 
tions present  on  a  carrier  flight  deck.   It  was  observed 
that  after  being  subjected  to  a  fire  which  resulted  in  the 
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epoxy  burning  away,  the  graphite  fibers  would  continue  to 
smolder  or  burn.   In  other  instances,  the  fibers  would  cool 
depending  on  the  flow  velocity  and  plate  thickness.   To 
understand  this  behavior  better,  a  mathematical  model  was 
formulated  to  predict  the  conditions  for  which  the  exo- 
thermic process  is  self-sustaining. 

The  air  flow  over  one  surface  produces  a  pressure 
differential  across  the  plate  which  induces  a  convective 
air  flow  through  the  porous  medium  governed  by  Darcy ' s  Law. 
As  a  result,  there  is  an  enhancement  of  internal  convection 
heat  transfer,  as  well  as  a  source  of  oxygen  for  combustion. 
The  heat  transfer  mechanisms  included  in  the  model  are 
(1)  conduction,  (2)  convection,  and  (3)  radiation.   In  addi- 
tion, non-volatile  combustion  is  included  in  the  energy 
balance  as  a  heat  generation  term  of  Arrhenius  type.   The 
oxygen  concentration  is  accounted  for  by  a  molecule-mass 
balance  which  includes  (1)  molecular  diffusion,  and  (2)  mass 
transfer  by  convection.   Consumption  of  oxygen  due  to  combus- 
tion is  accounted  for  by  a  term  similar  to  the  heat  generation 
term. 

The  formulation  of  a  one-dimensional  model  is  presented. 
All  properties  are  treated  as  temperature  dependent  in  the 
transient  analysis.   To  account  for  heat  transfer  mechanisms 
between  the  air  and  porous  graphite  medium,  the  air  and 
graphite  temperatures  are  treated  as  independent  variables. 
The  energy  balance  on  the  graphite  and  on  the  air,  and  the 
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molecule-mass  balance  on  oxygen  result  in  a  system  of  three 
nonlinear,  coupled  partial  differential  equations. 

Integration  of  the  field  equations  is  accomplished  by 
a  Galerkin  formulation  of  the  finite  element  method.   Due 
to  the  inherent  stiffness  of  the  field  equations  in  the  time 
domain,  a  modified  implicit-Gear  integration  scheme  developed 
by  Franke  [2]  is  used  for  a  more  efficient  solution. 
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II.   THEORY  AND  BACKGROUND 

A.   INTERACTION  OF  HEAT  TRANSFER  AND  COMBUSTION 

In  this  work,  the  combustion  model  proposed  by  Semenov 
was  adopted.   It  is  described  in  the  texts  of  Frank-Kamenetskii 
[3]  and  Vulis  [4] .   A  brief  discussion  of  those  features  of 
the  model  which  relate  to  the  present  investigation  will  be 
given.   Fundamental  to  the  model  is  the  relation  of  reaction 
rate  to  temperature,  and  the  interaction  of  the  heat  genera- 
tion and  the  heat  transfer  from  the  system.   In  general,  the 
reaction  rate,  r,  may  be  shown  by  the  Arrhenius  law  as, 


r   =   A(T,<j>)  exp(-E/RT) 

where  A   is  the  characteristic  time  of  the  chemical  reac- 
tion, E  is  the  activation  energy,  R  is  the  universal  gas 
constant,  T  is  absolute  temperature  and  <j>  is  the  concentra- 
tion of  oxygen.   The  concentration  appears  as  $   in  A(T,<j>) 
for  an  nth  order  reaction.   The  heat  generated  by  the  exo- 
thermic process  is  obtained  by  multiplying  the  reaction 
rate,  r,  by  the  enthalpy  of  formation  of  the  combustion 
products. 

In  Figure  1,  the  heat  generation,  q  ,  is  plotted  as  a 

y 

function  of  temperature.   The  curve  is  referred  to  as  the 
S-curve  for  apparent  reasons.   The  S-curves  are  distinguished 
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by  two  regions.   Lower  temperatures  and  reaction  rates 
characterize  region  I,  while  region  II  is  characterized  by 
higher  temperatures  and  reaction  rates.   In  region  I,  for 
low  reaction  rates,  there  is  an  abundant  supply  of  oxygen. 
The  reaction  is  controlled  by  the  temperature  and  the  region 
is  known  as  the  kinetic  regime.   Here,  the  reaction  rate 
increases  exponentially  with  increasing  temperature.   In 
region  II,  the  reaction  is  limited  by  oxygen  concentration 
and  is  known  as  the  diffusion  regime.   The  reaction's  weak 
dependence  on  temperature  in  region  II  is  shown  by  the  char- 
acteristic flattening  of  the  S-curve  in  Figure  1.   For  higher 
concentrations  of  oxygen,  the  diffusion  region  is  associated 
with  higher  reaction  rates.   The  interaction  of  heat  genera- 
tion, q  ,  and  heat  loss,  q„,  will  now  be  described, 
^g  ^£ 

In  addition  to  the  S-curve,  Figure  1  shows  two  heat 
loss  lines,  q  ,  and  qp?.   They  represent  the  heat  transfer 
by  convection  for  air  temperatures,  T-,  and  T2,  respectively. 
The  intersection  of  the  heat  generation  curve  and  the  heat 
loss  curve,  q   =  q  ,  defines  the  stable  quasi-stationary  point, 
A.   As  the  air  flow  temperature  increases  from  T-,  to  T2,  the 
quasi-stationary  temperature  increases  continuously  from  T, 
to  T_.   For  temperature  T  ,  q  '     is  tangent  to  q   and  an 
infinitesimal  increase  in  the  air  flow  temperature  results 
in  a  jump  of  the  reaction  temperature  from  TT  to  T_. .   At  point 
I,  which  is  defined  as  the  "critical  ignition  condition", 
the  reaction  moves  from  the  kinetic  to  the  diffusion  regime. 
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Temperature  TT  is  referred  to  as  the  ignition  temperature. 
Frank-Kamenetskii  states  the  ignition  temperature  of  graphite 
is  approximately  1100  to  1300  degrees  Celsius.   In  order  to 
keep  the  reaction  model  simple,  the  present  investigation  was 
limited  to  the  kinetic  regime.   Although  the  description  of 
the  heat  transfer  process  just  presented  is  for  convection 
only,  the  underlying  ideas  are  valid  for  the  conduction  and 
radiation  heat  transfer  mechanisms  as  well.   The  slope  of  the 
heat  loss  lines,  q  ,  is  equal  to  the  product  of  the  internal 
heat  transfer  coefficient  and  the  surface  area  of  fiber  per 
unit  volume.   The  slopes  would  not  appear  constant,  and  slight 
curvature  in  the  heat  loss  lines,  q  ,  is  apparent  when  radia- 
tion effects  are  included.   As  shown  in  Figure  1,  an  increase 
in  the  heat  transfer  rate  resulting,  for  example,  by  an 
increase  in  the  internal  flow  rate  would  result  in  a  lower 
graphite  temperature  for  a  given  air  temperature. 

One  reason  for  limiting  the  model  to  the  kinetic  regime 
was  due  to  the  jump  in  temperature  and  reaction  rate  at 
ignition.   Though  the  discontinuity  may  be  accounted  for,  the 
results  show  that  for  the  region  of  temperatures  approaching 
the  diffusion  regime,  the  heat  generation  will  dominate  the 
process.   The  system  temperatures  then  rise  rapidly.   Other 
factors  that  limit  the  model  to  the  lower  temperatures  are 
(1)  the  chemical  reaction  fixing  the  heat  of  combustion  (com- 
plex at  higher  temperatures)  and  (2)  the  rapid  consumption  of 
fibers  after  ignition. 
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B.   HEAT  TRANSFER  EQUATIONS  FOR  POROUS  MEDIA 

In  their  investigation  of  oil  recovery  from  underground 
reservoirs,  Green  and  Perry  [5]  developed  heat  transfer 
equations  for  the  solid  and  fluid  phases  in  a  porous  medium. 
Their  model  included  the  three  basic  mechanisms:   (1)  physical 
movement  of  the  fluid  which  carries  its  own  heat  capacity, 
(2)  conduction  of  heat  through  the  solid  and  fluid  phases, 
and  (3)  convective  heat  transfer  between  the  solid  and  fluid 
phases.   Radiation  heat  transfer  was  assumed  to  be  negligible, 
and  combustion  was  not  a  consideration.   The  differential 
equations  for  the  fluid  and  solid  phases 

2 
3Tf  3Tf       3  T- 

Fluid:   PfcfP—  =   -uPfcfPlF+  kfp— T-   hz(Tf-Ts)   (II. 1) 

o  X 

3T  32T 

Solid:   pscs(l-p)T^-  =   kg(l-p) ^-+hz(Tf-Ts)      (II. 2) 

3  x 

were  solved  numerically  by  the  finite  difference  method. 
In  equations  II. 1  and  II. 2,  c  is  the  specific  heat  at  constant 
pressure,  u  is  the  fluid  velocity  through  the  medium,  h  is 
the  coefficient  of  heat  transfer,  and  k  is  a  pseudothermal 
conductivity.   The  model  did  not  include  change  of  properties 
with  temperature.   Their  numerical  results  showed  good  agree- 
ment with  their  experimental  results.   Two  useful  conclusions 
derived  from  their  analysis  were:   (1)  both  conduction  and 
convection  heat  transfer  mechanisms  are  important  for  internal 
Reynolds  numbers  less  than  one,  and  (2)  for  values  of  the 
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1/2 

dimensionless  parameter,  £  =    (hz/k  )  /  (k_/p_cfu)  greater 

than  .342,  the  fluid  temperature  is  approximately  the  solid 
temperature.   Riaz  [6]  proposed  the  following  equation 
resulting  from  the  equal  fluid-solid  temperature  assumption, 


PSCs3t  +  pfCfU3^   =   kS  7-2  (II'3) 

a  X 


where  the  coefficients  are  defined  as  those  in  equations 
II. 1  and  II. 2.   For  the  present  model,  the  more  general 
approach  (i.e. ,  Green  and  Perry)  will  be  taken  in  the  formula- 
tion of  the  energy  equations . 


19 


III.   FORMULATION  OF  THE  FIELD  EQUATIONS 

A.   INTERNAL  FLOW  MODEL 

Referring  to  Figure  2,  the  porous  mat  may  be  represented 
by  a  flat  plate.   One  side  is  exposed  to  still  air,  and  the 
other  to  an  air  flow  of  velocity,  U  .   The  environment  on 
both  sides  is  at  ambient  temperature,  T  . 

The  first  step  in  developing  the  model  was  to  derive 
an  expression  for  the  flow  through  the  porous  plate.   Effects 
of  the  combustion  process  on  the  physical  properties  were 
neglected.   Specifically,  these  were  (1)  variations  in  the 
porosity  due  to  fiber  consumption,  and  (2)  effects  on  the 
viscosity  and  density  of  air  from  the  introduction  of  combus- 
tion by-products  into  the  flow.   Preliminary  calculations 
showed  that  the  porosity  and  density  of  the  material  did  not 
change  appreciably  until  temperatures  reached  approximately 
2000  degrees  Fahrenheit.   Further,  since  eighty  percent  of 
air  by  volume  is  essentially  inert,  it  was  reasonable  that 
its  properties  would  not  be  significantly  altered  during  the 
combustion  process.   Assuming  a  heterogeneous  material  of 
constant  porosity,  a  Reynold's  number  for  the  interior  flow 
is  defined  as 


Re.   =   p  u  d/y  (III.l) 

l     Ka  p 
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where  p   is  the  mass  density  of  air,  u   is  the  pore  velocity, 
d  is  the  diameter  of  the  fiber,  and  y  is  the  dynamic  viscosity 
Extensive  experimental  work  shows  that  Darcy's  law  for  flow 
through  a  porous  medium  is  valid  for  a  limited  range  of 
Reynolds  numbers.   Muskat  [7]  proposed  that  Darcy's  flow  is 
valid  for  Re.  <  1.   Scheidegger  [8]  points  out  that  a  number 
of  investigators  give  a  much  higher  upper  limit.   Preliminary 
calculations  of  Re.  for  the  model  showed  typical  magnitudes 
on  the  order  of  1.   Neglecting  the  gravity  term,  Darcy's  law 
for  flow  in  porous  media  is, 


m  dP  ,  _  _  _  ~  x 

u    = -3~  (III. 2) 

p        y  dx 


assuming  a  constant  pressure  gradient,  equation  III. 2 
becomes, 


m  AP  /ttt  r>\ 

u   =    =-  (III.  3) 

p       y  I 


where  AP  is  the  pressure  differential  across  the  plate 

thickness,  m  is  the  specific  permeability  of  the  medium, 

and  y  is  the  dynamic  viscosity.   In  equations  III.l  and  III. 3, 

p   and  y  were  treated  as  temperature  dependent, 
a 

Specific  permeability  is  a  measure  of  a  porous  medium's 
ability  to  allow  fluid  to  flow  through.   It  is  dependent  only 
on  geometry  and  the  physical  nature  of  the  medium.   It  has 
dimensions  of  length  squared.   Specific  permeability  is  often 
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called  the  hydraulic  conductivity  due  to  the  similarity  of 
Darcy ' s  law  to  that  of  Fourier's  law  for  heat  conduction. 
The  permeability  for  a  particular  medium  is  normally  measured 
by  experimentation.   However,  there  are  several  empirical 
models  that  may  be  used  to  obtain  values  for  permeability  that 
are  in  agreement  with  experimental  results.   For  the  idealized 
geometry  of  the  porous  medium  shown  in  Figure  2,  the  permea- 
bility was  obtained  by  a  serial  type  model,  equation  III. 4, 
proposed  by  Scheidegger  [8]  , 


m   =   -^p(6/x)2  (III. 4) 


where  the  porosity,  p,  is  defined  for  a  prous  material 
comprised  of  fibers  with  cylindrical  shape  as, 


p   =   1  -  j(d/D)2  (III. 5) 


and  D  is  the  thickness  of  a  ply.   Porosity  has  units  of 
volume  of  void  per  unit  volume  of  medium.   The  average  pore 
diameter,  5  was  defined  as, 

6   =   (s  +  D)/2 

for  the  geometry  in  Figure  2.   The  tortuosity,  x ,  is  a  measure 
of  length  of  travel  for  a  fluid  particle  per  unit  thickness 
of  the  medium.   Referring  to  the  geometry  of  Figure  2,  the 
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tortuosity  depends  on  the  ratio ,  d/D.   The  lower  limit 
occurs  when  d/D  equals  zero  and  for  the  upper  limit,  d/D 
equals  1.   Carman  [9]  presents  a  tab.le  of  tortuosities  for 
given  geometries   and  particle  shapes.   For  the  model, 
Carman's  suggested  value  of  t  =  1.4  was  used. 

To  obtain  the  pressure  differential,  AP,  across  the 
porous  plate,  Bernoulli's  equation  was  used.   The  following 
observations  showed  this  to  be  a  valid  assumption.   Schlichting 
[10]  points  out  that  for  the  ratio  of  u  /U   in  the  range  of 
.0001  to  .01,  the  effects  of  "blowing"  or  "suction"  on  the 
potential  flow  over  the  porous  plate  may  be  neglected.   A 
typical  value  of  u  /U   for  the  model  at  which  U  =  25  kt.  was 

■*  *  p '    00  00 

.0028.   For  steady  flow  over  a  flat  plate,  the  flow  field 
outside  the  boundary  layer  may  be  described  by  Bernoulli's 
equation.   This  is  a  direct  result  of  the  Navier-Stokes  equa- 
tion.  In  addition,  the  pressure  gradient  across  the  boundary 
layer  may  be  taken  equal  to  zero.   Therefore,  Bernoulli's 
equation  (eq.  III. 6)  may  be  used  to  obtain  the  pressure 
differential  across  the  plate, 

,  p       u 
—  dx  +  -pr-     =   constant.  (III.  6) 

pa       2 

The  parameters  in  equation  III. 6  are  defined  as  follows: 

P,  pressure;  p  ,  the  density  of  air;  U,  velocity  of  the  air 
a 

over  the  flat  plate.   For  the  model,  the  density  of  air  was 
approximated  by  the  Ideal  Gas  law  as 
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pa   =   P/RaTa  (III. 7) 


where  R   is  the  gas  constant  for  air,  T   is  the  absolute 
a  a 

temperature  of  air,  and  P  is  the  pressure.   Substituting 
this  into  equation  III. 6,  gives 


R  T        2 
/  (-^)dP  +  ^-  =   constant 


and  upon  integrating,  the  expression  becomes 


U2 
R  T  In  P  +  -=-  =   constant 
a  a        2 


or, 


2  2 

Ul  U2 

RaTal  ln  Pl  +  T"  '      RaTa2  ln  P2  +  T 


Letting 


P,  =  P  ,   U,  =  0,   T  ,  =  T  0  =  T  ,   P0  =  PT ,   U0  =  U 
1    <*>'   1    '   al    a2    •    2    l    2    < 


and  substituting  these  in  above,  yields 

°- 

R  T  In  P    =   R  T  ln  PT  +  -T- 
Rearranging  the  previous  expression  gives, 
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u2 

RaT„ in  «w   =  -f 


Taking  the  exponential  of  both  sides  and  solving  for  P 
results  in 


PT   =   P  exp  (-U2/2R  T  ) 
L       °°       °°    a  °° 


From  the  above  expression  and  noting  that  AP  =  PT  -  P^, 


AP  may  be  expressed  as 


AP   =   P  [exp(-U2/2R  T  )  -  1] 

oo      r      oo'     a   °° 


Darcy ' s  law  in  approximate  form  (eq.  III. 3)  for  the  model 
is 


m  AP 
P        M  L 


and  substituting  for  AP,  it  follows  that, 


P  5 

u    =   {E   "[i  -  exp(-U^/2R  T  ) ] }  (III. 8) 

poo       y  L       c        °°    a  °° 


where 


1  tw5\  2 
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The  pore  velocity,  u   of  equation  III. 8  is  for  an 
isothermal  medium  at  ambient  temperature,  T  .   Since  the 

oo 

model  is  to  investigate  the  transient  problem,  a  general 
expression  to  obtain  pore  velocity  at  different  temperatures 
had  to  be  developed. 

Consideration  of  the  continuity  equation  for  one-dimen- 
sional steady  flow, 


d(paup)/dx  =   0  (III. 9) 


becomes 


u  (dp  /dx)  +  p  (du  /dx)   =   0 
pa'        a   p 


or 


u  (dp  /dx)   =   -  p  (du  /dx) 
pa'  a   p 


Multiplying  through  both  sides  by  dx,  the  expression 
becomes, 


u  dp    =   -  p  du 
pa        a  p 


Separating  variables  and  integrating  both  sides,  gives 
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p 

u 

a 

/ 

dco 

f    P    ^a. 

00 

~~ 

p 

u 

roo 

poo 

and  the  integral  becomes 


ln(p  /p  )   =   ln(u  /u  ) 


Taking  the  exponential  of  both  sides,  gives 


p  /p    =   u  /u 
a  °°      p°°  p 


Rearranging,  the  pore  velocity  is  expressed  as 


u   =   p  u  /p  (III. 10) 


Referring  to  the  Ideal  Gas  law  for  the  density  of  air,  and 
substituting  this  into  equation  III. 10  gives 


u       (p  u   R  T  )/P  (III. 11) 

p       »  p°°  a  a  ' 


The  pore  velocity  is  now  a  function  of  air  temperature  and 
pressure.   However,  in  preliminary  calculations,  the  pressure 
difference,  AP  was  small  and  for  the  range  of  parameters  of 
interest,  it  did  not  exceed  five  percent  of  the  magnitude  of 
ambient  pressure.   Therefore,  as  a  simplification,  it  was 
assumed  that  P  equals  the  ambient  pressure  in  equation  III. 11, 
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and  the  expression  becomes 


u   =   (p  u  R  T  )/P 
p       °°  p°°  a  a    < 


Noting  that  p^  =  P^/R  Tw,  the  pore  velocity  may  be  expressed 
as 


u   =   u   T  /T  (III. 12) 

p      p°°  a   °° 


where  u    is  obtained  from  equation  III. 8.   T  and  T  are 
poo  ^  a      °° 

the  absolute  temperatures  of  the  air  at  a  point  in  the  interior 
of  the  medium  and  of  the  environment,  respectively. 

B.   ENERGY  EQUATIONS 

Here  we  consider  the  development  of  the  energy  equations. 
However,  before  continuing,  discussion  must  be  made  as  to 
the  approach  taken  for  energy  balances  in  porous  media.   In 
previous  work  by  Denbigh  and  Turner  [11] ,  and  by  Colladay 
and  Stepka  [12]  ,  it  was  suggested  that  the  temperature  of 
the  porous  solid  and  of  the  fluid  be  considered  equal.   This 
greatly  simplifies  the  formulation  of  the  problem.   However, 
under  certain  conditions,  as  shown  by  Green  and  Perry  [5], 
the  assumption  of  equal  temperatures  may  not  be  valid.   With 
recent  developments  in  numerical  techniques  permitting  the 
solution  of  non-linear  systems,  the  development  of  the  present 
model  for  heat  transfer  in  a  porous  medium  was  not  restricted 
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to  the  assumption  of  equal  temperatures.   Therefore,  the 
temperatures  of  the  air  and  of  the  porous  solid  are  treated 
as  independent  variables  in  the  present  formulation. 

The  model  is  based  on  two  assumptions.   First,  due  to 
the  small  thickness  of  the  graphite  fibers,  the  temperature 
across  each  individual  fiber  was  assumed  to  be  constant. 
Second,  for  a  similar  size  consideration,  the  temperature  of 
the  air  in  each  individual  pore  was  assumed  constant.   How- 
ever, the  temperature  from  fiber  to  fiber  and  from  pore  to 
pore  was  not  restricted,  and  was  allowed  to  vary  according 
to  the  heat  transfer  mechanisms  described  next. 

To  perform  energy  balances  on  both  the  porous  solid  and 
on  the  air,  a  differential  volume  of  porous  material  was 
segregated  into  respective  volumes  of  constituents,  that  is, 
graphite  fibers  and  air  (Figure  3) .   An  energy  balance  may 
be  done  on  each  volume  separately.   The  convention  used  for 
the  balance  is: 


Heat  into 
dV 


Heat 
Generation 


Heat  out 
of  dV 


Increase  in 
Internal 
Energy 


where 


dV   =   dx  dy  dz 
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1.   Fiber  Heat  Transfer  Equation 

Considering  the  differential  volume,  (l-p)dV  of 

graphite  fibers,  and  "smearing"  the  fibers  to  form  a  macro- 

scopically  homogeneous  material,  the  heat  transfer  mechanisms 

are  shown  for  the  x-direction  in  Figure  4.   The  heat  transfer 

mechanisms  are:   q   d,  heat  conduction  through  the  medium; 

q     ,  convection  from  the  fibers  to  the  air;  q   , ,  radiation 
^conv  ^rad' 

transfer  from  fiber  to  fiber;  AHr(T  ,  <f>)  heat  generation  rate 

y 

per  unit  volume.   AH  is  the  enthalpy  of  formation,  and 

r(T  ,<f>)  is  the  fiber  mass  consumption  rate.   This  term  will 
y 

be  discussed  in  detail  later. 

Representing  the  terms  in  Figure  4  by  Taylor  series 
expansions  and  combining  them  with  respect  to  the  convention 
previously  stated,  the  energy  balance  is 


cond 
V>nd  +   *rad  +   4Hr<Tg<*>  U-p)dxdydZ      =      qcond   +  — j-5—  dx 

3(2-r^  SI? 

+   q      a   +      .  dx   +   q  +   f£(l-p)dxdydz 

^rad  3x  ^conv        at        ^  J 


Subtracting  terms,  q    ,  and  q   ,  from  each  side  and 
*       '  ^cond     ^rad 

rearranging,  the  expression  becomes  (neglecting  the  higher 
order  terms) 


^cond  ,     ^rad 


dx ^^  -  q„     +  AHr(T  ,  <j>)  (l-p)dxdydz 

v  y 

(l-p)dxdydz  (III. 13) 


3x    *"     3x     Mconv        g 

3E 


3t 
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The  heat  transfer  terms  in  equation  III. 13  will  now  be 
considered. 

a.   Thermal  Conduction 

Fourier's  law  for  conduction  may  be  written  for 
the  incremental  volume  as 

3T 

q    ,   =   -  k  -T-2dydz  (III.  14) 

^cond        g  3x  * 

where  k   is  the  effective  conductivity  of  the  porous  solid, 
and  T   is  the  temperature  of  the  fibers.   There  are  a  variety 
of  models  for  the  effective  conductivities  of  porous  media. 
They  depend  largely  on  the  geometry  and  on  the  nature  of 
the  constituents.   For  the  cylindrical  shape  of  the  fibers, 
the  empirical  relation  proposed  by  Rohsenow  [13]  was  used 

k   =  k  U-e[l  +  Y  1  /9lnr~u~Y  } )]}   (III. 15) 

g      a         (1-Y)1/2        Y 

2 

for  y   <  1.   The  parameters  of  equation  III. 15  are  defined 

as  follows:   3  =  d/D  (refer  to  Figure  2),  y    is 


-   1r   A   i 

Y  "  ;cTT^aT] 


where  A  =  k  /k',  k   is  the  conductivity  of  the  air;  and  k1 
a.     g   a  J  g 

is  the  actual  conductivity  of  the  graphite  in  bulk  form. 
Since  k   depends  on  temperature,  k  was  also  treated  as 
temperature  dependent. 
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b.   Radiation  Heat  Transfer  in  Fibers 

Radiation  heat  transfer  in  the  model  was  repre- 
sented by  an  analog  to  Fourier's  law  of  conduction  as 

^rad  "  "  kr  T3?  ^  dz  (III. 16) 

where  k   is  an  equivalent  conductivity  of  the  fibers  due 
to  radiation.   Simple  assumptions  and  approximations  may  be 
made  to  show  this.   Assuming  air  to  be  transparent  to  radi- 
ation, and  treating  a  ply  of  graphite  fibers  as  an  infinite 
wall/  the  net  heat  flux  between  plies  may  be  written  as 


^ad  "  T2rfl(igx  -  *gx+dx>  (III. 17) 


where  T   and  T  x+dx  are  the  respective  absolute  ply  tempera- 
tures, e  is  the  emissivity  of  the  graphite,  and  a  is  the 

Stefan-Boltzman  constant.   Treating  T   as  the  variable 

gx 

temperature,  q"  ,  may  be  expanded  in  a  Taylor  series  about 
T      .   Neglecting  the  higher  order  terms,  the  series 
expansion  may  be  written  as 


a«  =   ^-(T4      -  T4     )  +  i°-l  T3     (T   -  T      ) 

Mrad      2-ev  gx+dx   gx+dx'    2-e   gx+dxv  gx    gx+dx; 


Simplifying,  the  above  equation  becomes, 


„     =   4o_e_  p  (£  -i  )  (III. 18) 

^rad      2-e   gx+dx   gx    gx+dx 


32 


Taking  Fourier's  law  for  heat  transfer, 

dT 

q"  *    =    -  * 


rad        r  dx 


dT 
and  approximating  -3-2.  with 


dT2     ^    Tgx+dx  "  Tx 


dx      Ax  Ax 


["  ,  becomes 


(T   _   -  T   ) 
_     _JE+dx gx_ 

^rad        r       Ax 


Multiplying  the  right  side  of  equation  III. 18  by  Ax/Ax, 
equation  III. 18  and  equation  III. 19  may  be  equated 


(T      -  T  )  (T    —  -  T 

gx+dx    gx   _     4aeAx  ,13       gx+dx 


^rad        r      Ax  2-s    gx+dx      Ax 


Cancelling  terms,  k  becomes 


gx 


k   =   ioeAx  T3  (in. 20) 

r       2-e    gx+ax 


where  Ax  is  now  equal  to  5,  the  average  pore  diameter. 
From  the  close  spacing  of  the  fiber  layers,  the  temperature 
difference  will  be  small  as  compared  to  the  magnitude  of  the 
temperature.   Noting  this  the  average  absolute  temperature  of 
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the  fibers.  T  may  be  substituted  for  T   , ,  .   The  equivalent 
g    J  gx+dx 

radiation  conductivity  expression  becomes 


k    =   4°Si  T3  (III. 21) 

r      2-e   g 


and  the  radiation  heat  transfer  from  fiber  to  fiber  may  be 
represented  by, 


•rad  "   ""r^f  =  "  TT  Tg  S  ^   dz     (III'22) 


If  the  temperature  gradient  in  the  fibers  is  not 

large,  then  equation  III. 22  will  be  a  good  approximation  of 

the  radiation  heat  transfer  between  the  fibers. 

c.   Convection  Heat  Transfer 

The  convection  heat  transfer  term  in  Figure  4 

was  treated  in  a  similar  manner  to  that  of  the  one-dimensional 

fin  equation.   q     was  introduced  as 
^         ^conv 


q       =   h.  (T  -  T  )  dA  (III. 23) 

nconv      l   g   a 


where  T   is  the  temperature  of  the  graphite,  T   is  the 
g  a. 

temperature  of  the  air,  dA  is  the  surface  area  of  graphite 
in  the  differential  volume,  and  h.  is  the  internal  convection 
heat  transfer  coefficient.   An  empirical  expression  in  the 
form  of  a  Colburn  j -factor  was  developed  by  Yoshida, 
Ramaswami,  Hougen  [14].   This  was  used  to  determine  h., 
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h  •    c  V     9  /->  _  n 

J   =   — i_[_£_]^/J   =   .91  Re~'0±V  (III. 24) 

C,  G    K  O 

bo   a 


for  Re  <  50.   c,  is  the  effective  specific  heat  of  the  porous 


medium  and  is  defined  as, 


cb  =  pca  +  (1-P)cg 


where  c   is  the  specific  heat  of  air  at  constant  pressure; 
a 

c   is  the  specific  heat  of  graphite;  and  p  is  the  porosity 
y 

as  defined  by  equation  III. 5.   G   is  a  pseudo  mass  velocity 

defined  as  G   =  p  u  p.   Also  in  equation  III.  24,  \i    is  the 
o    a  pr  ^ 

viscosity  of  air,  and  k   is  the  conductivity  of  air.   The 
Reynolds  number  appearing  in  equation  III. 24  is  defined  as 


Re  =   G  /z  u    <y 
o      o 


and  is  not  the  same  as  Re.  defined  previously  for  Darcy ' s 
law.   z  is  the  surface  area  of  graphite  fibers  per  unit 
volume,  and  is  determined  from  geometrical  considerations 
which  will  be  described  shortly.   The  parameter,  ¥  ,  is  a 
dimensionless  shape  factor  which  depends  on  the  geometry  of 
the  fibers.   For  a  cylindrical  fiber  shape,  ¥   is  equal  to 
.91.   From  the  assumption  of  constant  porosity,  z  remains 
constant.   However,  the  air  properties  are  temperature  dependent; 
therefore,  h.  is  allowed  to  vary  with  temperature. 
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d.   Heat  Generation  Rate 

The  reaction  rate,  r(T  ,<f>)/  whether  based  upon 
the  Collision  Theory  or  on  the  Theory  of  Absolute  Reaction 
Rates,  results  in  an  expression  of  Arrhenius  type, 


r   =   A(T,<}>)  exp(-E/RT) 

where  E  is  the  activation  energy,  R  is  the  universal  gas 
constant,  and  T  is  absolute  temperature.   Although  E  is  well 
defined  for  the  combustion  of  graphite,  A(T,<j>)  is  a  function 
of  temperature  and  concentration  and  depends  upon  chemical 
kinetic  theory.   Development  of  a  general  model  for  the  chemi- 
cal kinetics  is  beyond  the  scope  of  this  work.   A  less  ele- 
gant, but  useful  approach  of  using  a  relation  obtained 
experimentally  was  taken.   A  relation  derived  from  the  work 
of  Parker  and  Hottel  [15]  was  used  to  predict  the  combustion 
rate  of  graphite  fibers, 

a    n-u   m7  „   ml/2         ,-39,883*     lbm     ,___  oc, 
r   =   4.014  x  10   R  T  '       y  <$>    exp  ( r-* )   — ^ (III. 25) 

A      g  T       ft  -hr 

g 

The  development  of  equation  III. 25  from  the  original  work 
may  be  found  in  the  appendix.   This  expression  assumes  a 
simple  first-order  reaction  for  the  combustion  of  graphite 
in  air.   However,  Frank-Kamenetskii  [3]  has  further  refined 
Parker  and  Hottel' s  expression  to  yield  for  the  present 
model, 
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a    at  a      ir>7   „      ml/2         ^n            ,-39,883,           Ibm  /-n-T    nc\ 

r      =      4.014  x  10      R      T  '       y    <j>      exp  ( '- )       — = (III. 26) 

A  g  T      ft  -hr 

g 

The  combustion  is  now  characterized  as  an  nth  order  react- 
tion,  where  n  is  in  the  range  of  1/3  to  2/3.   This  form  of 
the  reaction  rate  expression  better  approximates  Parker  and 
Hottel's  experimental  data.   For  the  analysis,  n  was  chosen 
to  be  1/2.   However,  it  should  be  pointed  out  that  the  behavior 
of  the  reaction  varies  significantly  within  the  range  of  n 
proposed  by  Frank-Kamenetskiii. 

In  limiting  the  model  to  the  kinetic  regime  of 
combustion,  it  was  assumed  that  the  dominant  chemical  reac- 
tion is 


C  +  o2   ■  co2 


This  is  an  idealized  treatment  of  the  reaction  that  actually 
takes  place.   Kanury  [16]  states  that,  the  actual  reaction 
is  more  complex,  yielding  various  concentrations  of  carbon 
dioxide  and  carbon  monoxide,  depending  on  the  temperature. 
Noting  the  abundant  oxygen  supply  present  in  the  kinetic 
regime,  it  is  reasonable  that  the  leaner  oxide  is  the  preva- 
lent by-product  of  combustion.   The  fuel-oxygen  ratio,  f, 
for  this  mechanism  is  12/32.   Carbon  monoxide  may  be  formed 
when  there  is  little  free-oxygen  present,  such  as,  during 
combustion  in  the  diffusion  regime.   The  stoichiometric 
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fuel-oxygen  ratio,  f  =  12/16,  for  this  mechanism  is  naturally 
greater.   To  obtain  the  heat  generation  rate,  the  reaction 
rate,  r,  was  multiplied  by  the  enthalpy  of  formation,  AH, 
of  carbon  dioxide.   AH  has  units  of  Btu  per  lbm  of  graphite 
consumed. 

e.   Change  in  Internal  Energy 

The  last  term  of  the  graphite  energy  equation, 

3E 
the  rate  of  change  of  internal  energy,  —,  results  by 

setting  E  =  p  c  T  .   Thus, 

g  g  g 

_  3T 

|£(l-p)dxdydz   =   (l-p)p  c  -jg-   dxdydz  (III. 27) 

where  p   is  the  density  of  graphite  in  bulk;  c   is  the 
specific  heat  of  graphite;  and  p  is  the  porosity. 

The  heat  transfer  mechanisms  for  the  graphite 
fibers,  equations  III. 14,  III. 16,  III. 23,  III. 26,  III. 27  are 
substituted  into  equation  III. 13  which  upon  dividing  through 
by  dV  =  dxdydz  yields  the  heat  transfer  equation  for  the 
fibers, 


(1-p'  k%  t£>  +  (1-p'k<kr  S>  "  hi  af<VV 


3T 


+  (l-p)AHr(T  ,<$>)      =   (l-p)p  c 


J2 


g'r/     v~  r  g  g  3t 

The  coefficient,  (1-p) ,  accounts  for  the  fact  that  not  all 
of  the  volume  of  porous  medium  is  comprised  of  graphite. 


38 


Dividing  through  by  (1-p) ,  combining  the  first  two  terms, 
and  defining  dA/dV  as  z,  the  surface  area  of  graphite  per 
unit  volume,  the  expression  becomes 

3T      h.z 

4-[  (k  +  k  )-t-2-]  -  /11  t  (T  -  T  )  +  AHr(T  ,<j>) 
3x   g   r  3x     (1-p)   g    a         g  T 

3T 

=   P  C   -rq?  (III. 28) 

g  g  3t 

The  parameter,  z,  is  obtained  from  geometry  considerations. 
For  a  cylindrical  fiber  shape,  the  expression  for  z  is, 


z   =   7rd/2D2  (III. 29) 


The  one-half  factor  in  equation  III. 29  accounts  for  the  un- 
certainty of  the  actual  amount  of  exposed  fiber  surface  area, 
which  resulted  from  the  initial  combustion  of  epoxy.   The 
non-dimensionalization  of  equation  III. 28  is  presented  in 
Appendix  A. 

2.   Internal  Flow  Heat  Transfer  Equation 

The  energy  balance  on  the  internal  flow  is  obtained 
by  a  procedure  similar  to  the  energy  balance  on  graphite. 
Figure  5  presents  a  "smeared"  differential  volume  of  air, 
pdV,  which  shows  the  heat  transfer  mechanisms.  The  heat 
transfer  mechanisms  are:   q    ,,  heat  conduction  through  the 

air;  q     ,  convection  from  the  fibers  to  the  air;  m  h, 
xonv  a 
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energy  transport  due  to  the  air  flow,  where  h  is  the 
enthalpy  of  the  air . 

Performing  the  energy  balance  (neglecting  the  higher 
order  terms)  yields, 


^cond 


q    ,+mh  +  q       =   q    ,  +  r dx 

^cond    a    ^conv     ^cond     3x 


+  m  h  +  m  |$  dx  +  If  pdxdydz  (III. 30) 

a     a  o  x      a  x. 


Cancelling  terms  and  rearranging,  the  expression  becomes, 


cond  ,     •   3h  ,  3E,,j       ,TTT  „, 

-IF"  dx  "  ma  3x"  dx  +  qconv   =   3t  PdxdYdz      (HI.  31) 


a.   Thermal  Conduction 

As  before,  q    ,  was  replaced  by  Fourier's  law 
of  heat  conduction  in  the  form, 


>=ond  "  "  ka  ^T  d*dz  (III-32) 


where  k   is  the  conductivity  of  air. 
a 

b.   Energy  Transport  by  Flow 

Using  the  perfect  gas  assumption  for  air,  the 
enthalpy  h  was  expressed  as, 


dh  =   c   dT 
a   a 
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where  c  is  the  specific  heat  of  air  at  constant  pressure. 
The  specific  heat  is  treated  as  constant  since  it  does  not 
vary  appreciably  over  the  range  of  temperatures  associated 
with  the  problem  under  investigation.   It  follows  that, 


,;      3T 

~   =   C   -r-^  (III. 33) 

3x      a  3x 


Assuming  the  effects  of  the  combustion  by-products  on  air 
are  negligible,  continuity  considerations  give  the  mass 
flow  as, 


m   =   p   u  dydz  (III. 34) 

a      a   p 

OO   *■  00 


where  the  density  of  air  and  the  pore  velocity  are  evaluated 

at  ambient  conditions . 

c.   Convection  Heat  Transfer 

The  convection  heat  transfer,  q     given  by 

nconv  3      2 

equation  III. 23,  with  h.  obtained  from  equation  III. 24, 
becomes 


^conv  =  VVTa)dA  (II1-35) 


d.   Change  in  Internal  Energy 

For  a  perfect  gas,  the  internal  energy  is, 


dE  =   p  c  dT 
a  a   a 
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where  c   is  the  specific  heat  of  air  at  constant  pressure. 
From  this  it  follows  that 


3E  3Ta 

r-r-  p  dxdydz   =   p  c   >— —  p  dxdydz        (III.  36) 

()t  a  a  a x 


Substituting  equations  III . 32-111 . 36  into  equation 
III. 31,  the  heat  transfer  equation  for  the  internal  flow 
becomes, 


3T  3T 

"&(ka  -alf  >  P<*xdydz  -  maca  -^  p  dxdydz 


3T 

+  h.  dA(T  -T  )   =   pc,  -r~  P  dxdydz 
i      g   a        a  o  u 


As  in  the  energy  equation  for  the  fibers,  the 
porosity,  p,  accounts  for  the  fact  that  some  of  the  volume  of 
the  porous  medium  is  not  air.   Dividing  the  last  equation 
through  by  pdxdydz  and  letting  z  denote  dA/ dxdydz ,  yields 


3T  9T     h.z 

4~(k   -r-2.)  -  m  c   -r#  +  -^-(T  -T  ) 
3x  a  3x      a  p  3x     P   g   3. 


m     paca  TF  (III. 37) 


Non-dimensionalization  of  equation  III. 37  may  be  found  in 
Appendix  A. 


42 


3.   Oxygen  Transport  Equation 

The  final  consideration  in  the  formulation  of  the 
model  is  the  mass  transport  of  oxygen.   The  mass  transport 
equation  is  obtained  by  a  mass  balance  on  a  differential 
volume  of  porous  medium.   In  accordance  with  convention, 
the  mass  balance  is  given  by 


Oxygen  into   _   Oxygen  out      Oxygen        Oxygen 

dV  of  dV     Consumption   Accumulation 


Figure  6  presents  a  differential  volume,  pdV,  showing  the 
relevant  molecule-mass  transport  mechanisms.   The  mass 
transport  mechanisms  are:   m_  is  the  molecular  diffusion; 
m  is  the  convective  transport  due  to  mass  flow;  ^-r(T  ,  <f>) 
is  the  consumption  of  oxygen  due  to  combustion;  r  is  the 
reaction  rate,  and  f  is  the  stoichiometric  ratio  of  the 
reaction.   Representing  the  terms  in  Figure  6  by  Taylor 
series  expansions,  the  molecule-mass  balance  becomes  (neg- 
lecting higher  order  terms) , 


3mB  3™c 

mB  +  ™c  "  ^B  +  ax  dx  +  ™c  +  TT  dx 


+  ir(T  ,<|>)pdxdydz  +  |£pdxdydz  (III.  38) 


Upon  cancelling  terms  and  rearranging,  equation  III. 38 
becomes 
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3  m  3m,  a  . 

-^-=-  dx   -   -r-£  dx   -   ±r(T    ,  «>)pdxdydz      =      |A  pdxdydz       (III.  39) 

o  X  dX  X     g  " "~ 


where  <j>  is  the  concentration  of  oxygen.   The  individual  mass 
transfer  terms  in  equation  III. 39  will  now  be  developed, 
a.   Molecular  Diffusion 

Molecular  diffusion  nu,  is  given  by  Fick's  law, 


mB   =   -  B  ||  dy  dz  (III. 40) 


where  B  is  the  diffusivity  of  oxygen  into  the  porous  medium. 

The  diffusion  of  gases  in  porous  media  is  limited 
by  two  factors.   One  is  the  collision  of  the  gas  molecules 
with  other  molecules,  and  the  second  is  the  collision  of  the 
gas  molecules  with  the  solid  pore  walls.   This  second  phenom- 
enon is  discussed  by  Bennett  and  Myers  [17]  and  is  referred 
to  as  Knudsen  diffusion.   To  establish  which  behavior  dom- 
inates, the  mean  free  path  of  the  gas  (in  this  case,  oxygen) , 
must  be  calculated.   Using  the  expression  given  by  Treybal 
[18], 


w   =   (3.2y/P) (RT/2ug  M)  1/2 


where  g   is  the  gravitational  constant,  the  value  obtained 
for  mean  free  path,  w,  is  smaller  than  0.16.   Thus,  the 
collision  of  the  oxygen  molecules  with  those  of  air  will 
restrict  the  diffusion  process.   On  the  contrary,  if  the  pore 
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diameter  is  smaller  than  the  mean  free  path,  then  the  colli- 
sion of  the  molecules  with  the  solid  walls  will  be  the 
limiting  factor.   In  the  intermediate  case,  the  mean  free 
path  and  the  pore  diameter  are  of  the  same  order  of  magni- 
tude.  In  this  case,  Bennett  and  Myers  [17]  proposed  the 
following  expression  be  used  to  obtain  a  value  of  diffusivity, 

r  ■   w-  +  1r  (III-41> 

R       g     w 

where  B   is  the  resulting  diffusivity;  B   is  the  diffusivity 

based  on  molecule-molecule  collisions;  and  B   is  the  diffusivity 

w  ■* 

based  on  molecule-pore  wall  collisions. 

For  the  problem  under  investigation  here,  prelimi- 
nary calculations  of  the  mean  free  path  of  the  oxygen  showed 
that  diffusion  is  primarily  of  the  inter-molecular  collision 
type.   With  respect  to  this,  and  knowing  that  molecular 
diffusion  is  highly  temperature  dependent,  the  diffusivity, 
B*  was  obtained  from  an  expression  presented  by  Gilliland 
[19]  . 

J.3/2 
B*   =   435.7  -a—(V^/3+vJ/3)"2(M"1+MQ^)1/2   cm2/s 

(III. 42) 

where  P  is  the  absolute  pressure;  V  and  Vno  are  the  molecular 

a      u  <* 

volumes  of  air  and  oxygen  respectively;  Ma  and  Mno  are  the 

a      u  z 

molecular  weights  of  air  and  oxygen,  respectively,  and  B* 

is  the  diffusivity.   For  this  expression,  T   is  in  degrees 

a 
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2 

Kelvin,  and  B*  has  the  units  of  cm  /s.   As  a  result  of 

tortuosity  and  porosity  affecting  the  process,  Denbigh  and 
Turner  [11]  proposed  that  the  effective  molecular  diffusivity 
for  a  porous  medium  is  given  by 


Be   =   pB*/x. 


The  porosity  p  accounts  for  the  void  in  a  differential  cross- 
section. 

b.   Convective  Transport 

The  convective  term  in  equation  III. 39  is  given 
by 


m   =   pu  4>dydz  (III. 43; 


The  pore  velocity  u   is  a  temperature  dependent  variable, 
and  is  so  treated. 

c.   Consumption  Rate 

The  reaction  rate  term,  r  (T  ,  <j>)  was  discussed 
previously  as  equation  III. 26.   The  consumption  rate  of 
oxygen  is  determined  by  multiplying  r  by  1/f. 

Substituting  equations  III. 40  and  III. 43  into 
equation  III. 39,  the  oxygen  transport  equation  becomes, 


t — ( —  -r±)  pdxdydz  -  - — (u  <j>)  pdxdydz 
3x  t   3x   r    ■*      3x   p    * 


1  3  d> 

(f)r(T  ,  <J>)  pdxdydz   =   |i  pdxdydz 
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Dividing  both  sides  by  pdxdydz,  and  letting  B  =  B*/x,  the 
expression  becomes 


k<*&  -IsV  "  <?>r<V*>    -   H-        (III-44) 


Expanding  the  second  term  of  equation  III. 44, 


(u  di )   =   u   — —  +  £- 


and  noting  that  u  depends  on  T   from  equation  III. 12,  the 
expression  yields 


fx'V'   =  %£  +  ^TZ*   •  «"•«' 


3u  /3T   is  obtained  from  equation  III. 12. 
Pa  M 

Substituting  equation  III. 4 5  into  equation  III. 44, 
the  oxygen  transport  equation  becomes, 

k*&    --pH-S^T*  -  <£>r<V*>   -||     (III. 46) 


In  contrast  with  the  heat  transfer  equations,  mass  balance 
does  not  depend  on  porosity.   However,  porosity  does  influence 
the  oxygen  concentration  by  appearing  in  the  boundary  condi- 
tions.  Non-dimensionalization  of  equation  III. 46  is  pre- 
sented in  Appendix  A. 
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IV.   BOUNDARY  CONDITIONS 

It  was  intended  that  boundary  conditions  be  imposed 
which  simulate  the  composite  plate  "burn"  experiments  of 
Fontenot  [1] .   However,  for  a  one-dimensional  model,  the 
boundary  conditions  can  only  be  approximated  at  best.   The 
difficulty  is  associated  with  the  existence  of  momentum  and 
thermal  boundary  layers  on  the  exterior  surface  of  the  por- 
ous plate.   Complications  also  arise  from  the  "blowing  and 
suction"  effects  which  are  caused  by  the  flow  through  the 
plate. 

The  following  boundary  conditions  provide  a  reasonable 
model  for  the  one-dimensional  problem, 


3  T 

(kg  +  V  l£  =   hl(Tg~TJ  +  a£(ig  "  *t       at   x  =  0   (IV. 1) 

9T  ~4  -4 

(kg  +  V    Tf  =  "W^    "   a£(g  "      -)      at     x  =  L      (IV'2) 


T=T  atx=0  (IV. 3) 

a  °° 


3T  3T 

3X  3X 


at      x   =   L  (IV. 4) 


B  14     ■      u^(*    ~  P*J         at      x   =    0  (IV. 5) 

a  j\  p  °° 


|^=0  atx=L  (IV. 6) 

a  X 
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Boundary  conditions  IV. 1  and  IV. 2  account  for  convection 
and  radiation  heat  transfer  from  the  porous  solid.  With 
an  influx  of  air  at  x  =  0 ,  the  heat  transfer  coefficient, 
h, ,  will  depend  on  the  magnitude  of  the  flow  entering  the 
plate.  Neglecting  "blowing  and  suction"  effects,  several 
empirical  expressions  based  on  free  convection  were  used  to 

obtain  a  heat  transfer  coefficient  at  x  =  0.   These  expressions 

2 

yielded  a  range  of  values  for  h,  from  1  to  3  Btu/hr-ft  «F. 

As  an  approximation,  h,  was  taken  as  the  average  value,  2.0 

2 

Btu/hr«ft  *F.   However,  choosing  h..  as  any  value  in  the 

above  range  did  not  affect  the  results  of  the  analysis. 

The  magnitude  of  the  forced  convection  heat  transfer 
coefficient  at  x  =  L,  h2,  varies  in  a  direction  parallel  to 
that  of  the  external  flow,  Ua.   In  addition,  h2  depends  on 
the  efflux  of  the  gas  at  the  surface.   To  simplify  the  analy- 
sis, h2  was  approximated  by  the  relation  for  a  smooth  flat 
plate  given  by  Holman  [20]  as 


h2   =   .664  ^4  Pr1/3  Re1/2  (IV. 7) 


for  laminar  flow,  where  Pr  is  the  Prandtl  number  and  L*  is 

a  reference  length,  arbitrarily  taken  as  unity.   For  turbulent 

flow, 


h2      =      £§■  Pr1/3(.037   Re*8    -    850)  .  (IV. 8) 
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Kays  [21]  provides  an  alternate  scheme  for  treating  boundary 
layers  and  "blowing  and  suction"  effects. 

Boundary  condition  IV. 3  accounts  for  the  air  entering 
the  porous  plate  at  ambient  temperature.   The  Dankwerts 
boundary  condition  proposed  by  Riaz  [6]  which  accounts  for 
conduction-convection  interaction,  is  a  more  reasonable 
simulation. 

The  boundary  condition  at  x  =  L  for  the  air  heat  trans- 
fer equation  appears  to  be  a  reasonable  one.   It  follows  the 
behavior  observed  in  the  interior  of  the  plate  and  does  not 
fix  the  temperature  of  the  air  at  the  wall.   The  boundary 

condition,  3T  /3x  =  0,  used  initially  provided  results  simi- 
a 

lar  to  those  obtained  using  equation  IV. 4. 

Boundary  conditions  IV. 5  and  IV. 6  are  Dankwert  conditions 
[22]  for  flow  through  porous  media.'   A  convincing  discussion 
for  the  Dankwerts  conditions  is  given  by  Bischoff  [23] .   A 
brief  summary  of  the  discussion  is  presented  in  Appendix  D. 
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V.   FINITE  ELEMENT  METHOD 

A.   GALERKIN  FORMULATION 

A  Galerkin  formulation  of  the  Finite  Element  Method  was 
used  to  obtain  solutions  of  the  field  equations.   A  conven- 
ient form  of  the  field  equations  was  used  in  the  formulation 
These  equations  are  presented  in  Appendix  A. 
graphite  energy  equation 


3  0 
k(Al  lf»  "  X2(S-r)  +  A3r(Tg'*)   "   X4  "    (V-1) 


air  energy  equation 


k^i*^  -^+  Ve-r>    -   u4^    (v-2> 


oxygen  diffusion  equation 


t—  (w,  t— )  -  co~ u-  t—  *-w.r(T  ,A)   =   go-  — -       (V.3) 

3n   1  3n    2  3n    3  3n     4   g'r       5  3t 


where  the  coefficients  are  defined  in  Appendix  A. 

Field  equations  III. 28,  III. 37  and  III. 46  subject  to 
boundary  conditions,  IV. 1  -  IV. 6,  and  initial  conditions 
define  the  problem.   The  closed  domain  (0,L)  was  partitioned 
into  (n-1)  contiguous  elements  of  variable  length  l.,    i  =  1, 
. . . , (n-1) .   This  defines  an  n  nodal  point  model.   The  three 
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non-dimensional  field  variables  0,  r,  $  were  approximated 
by 

e(n,t)   =  ^1(n,t)      =     lGLMQ±(t)  (V.4) 

r(n,t)  =  4'2(n,t)  =  EG^TOi^Ct)      (v. 5) 

*(n,t)   =  ,?3(n,t)   =   [Gi(n)$i(t)      (V.6) 


where  G.,  for  i=l,  ...,  n,  is  a  set  of  specified  basis 
functions  with  local  support,  and  the  sets  { 0 . , r . , $ , ; i=  1,  ..., 

J»    -!•    JL 

n}  are  the  solution  coefficients  to  be  determined.   The  G. 

were  selected  to  satisfy  the  condition  G. (n-)  =  5. .,  where 

i  l   j      ID 

the  kronecker  delta  6 . .  is  defined  by  6 . .  =  1  for  i  =  j , 
and  6. .  =  0  for  i  ^  j.  As  a  result,  0,  r,  and  $  are  the 
values  v.  ,  v.,    v      at  the  nodal  points  (i.e.,  0.(t)  =  li'1(n-/t)). 

Linear  interpolation  functions  (Figure  7)  were  used 
as  the  basis  functions.   These  are  the  lowest  polynomial 
functions  which  provide  the  necessary  function  continuity. 

As  a  measure  of  error,  a  residual  function,  r.,  is 
defined  for  each  field  equation  by 


r±(n,t)   =   A±m  -  V±   i  =  1,  2,  3      (V.7) 


where  A.  denotes  the  spatial  operator  of  the  ith  equation. 
For  convenience,  the  following  convention  for  differentiation 
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is  adopted, 

3(  ) 


3n 


=   (  )  ' 


32(  )   =   ,  Y„ 

3n 


at     v  ; 


For  field  equations  V.l,  V.2,  and  V.3  the  residuals  are 


n  n 

r,     =     U,(e,r)     I     G.'e.]'   -  x9(r)      \     G.  (0.   -  r . ) 
1  1    -    -       >1      1   1  2.^11  1 


n 
+      A3r(0,$)    -  X4      J      G.^  CV.8) 

i=l 


■2  -  [vi(e>  .1  w  -  -2  .z.  Giri +  *3{v  .!.  Gi(ei-ri' 

1=1  1=1  1=1 


n 

4' 


„(r)      y     g.t-  (v.9) 

.  u  1  11 


i=l 


r3     =      [«    (r)      I     G!*.]'    -   «2(r)      I     G!$.-o>    (r)      V     G.r.G.*. 

i=l  i=l  i=l 


n 

u>,r(0,$)    -   u.      J      G.$.  (V.10) 

4      "    -  5   i=l      x   1 
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where  the  coefficients  multiplying  the  response  variables 
are  themselves  functions  of  the  response  variables,  and  thus 
the  equations  are  nonlinear.   In  accordance  with  the  Galerkin 
method,  the  final  system  of  ordinary  differential  equations 
was  obtained  by  setting  each  residual,  r.,  orthogonal  to 
each  basis  function,  G.,  that  is 

1 
/   G.r.dn   =   0   i  =  l,2,...,n;         (V.ll) 
0   X  J  TO 

j    =   1,2, ... ,n 

The  3n  ordinary  differential  equations  given  by  equations 
V.ll  retain  the  character  of  the  original  set  of  partial 
differential  equations.   Thus  linear  field  equations  trans- 
form to  matrix  operators  and  nonlinear,  coupled  field  operators 
become  nonlinear,  coupled  algebraic  operators.   Incorporation 
of  the  boundary  conditions  resulted  in  3n  nonlinear,  coupled 
ordinary  differential  equations 

A(t)  <F(0,r,$)  +  F(t)   =   B(t)  I  (V.12) 

subject  to  initial  conditions,  where  B  is  a  3n  x  3n  matrix, 
A  is  the  operator  associated  with  the  field  operator  A .  in 
expression  V.7,  and  F(t)  is  an  excitation  vector. 
Adopting  the  convention, 


n 

<Gi>  {».}   =    I      G.*. 

i       i=1   i  i 
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and  applying  the  operation  of  expression  V.ll  with  an  inte- 
gration by  parts  on  the  second  order  derivatives  gives, 


(G.  }A,<G.  >'  {0} 
ill 


1  1 

-    X,       /    {G.  }  '<G.>  'dnO} 
0  1       0      x  3 


1  1 

-    A9      /    {G. }<G.>dn(0}    +    Xj      S    {G. }<G.>dn(r}  (V.13) 


1  1 

+   x3      I    ^Gi>dnr (0,$)      =      A4      /    {G.^  <G . >dn {0 } 


(G, }v. <G.> ■ {r} 


ll 

-   vx      /    {Gi},<G.>'dn{r} 


1  1 

v2   /  {Gi><G.> 'dn{r}  +  v3   |  {G. }<G.>dn(0}       (V.14) 


0 


0 


1  1 

v   /  {G.}<G. >dn(r}  =  v4  /  {G. }<G.>dv{r} 


o 


o 
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{G.  }u>.  <G.>  '  {$} 
l      1      J 


11 

-  o>.       /    {G. } ,<G.>,dn{$} 
1       ;  n       l  1 


-  u -      f    {G. }<G.> 'dnU}    -    u,      /    (G. }<G.> ' {r }<G . >dn{* } 


1  1 

-  a).   /  {G.  }dnr  (0,$)   =   wc   /  {G .  }  <G .  >dn  { $}       (V.15) 


The  first  term  in  each  of  the  above  expressions  is  a  boundary 
term  which  permits  the  incorporation  of  natural  boundary 
conditions  which  will  be  shown  in  Section  V.B.   Coefficients 
X,  v,  u  are  variable-dependent  properties,  and  were  taken 
as  the  average  value  of  the  properties  over  an  element.   In 
the  limit,  as  the  elements  get  smaller  (i.e.,  n  -*■   «) ,  the 
average  values  converge  to  the  exact  values . 

Inspection  of  expression  V.13,  V.14,  V.15  shows  the  five 
operators, 


1 

/  {G.}'<G.>*dn  (V.16) 

«   i    i 


1 
/  {G. }<G.> 'dn  (V.17) 
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1 

/  {G.}<G.>dn  (V.18) 


/  {G.  }<G.>,{lF}<G.>dn  (V.19) 

0   1    ^        ^ 


1 
/  {Gi)dn  (V.20) 


To  formulate  these  operators,  the  global  shape  function,  G., 
was  defined  on  the  local  level  by 


Gi  =  g{i"1)  0  g^i}  (v.2i) 


where  g,  and  g_  were  defined  by 


(1  -  j-)      for  £  in  element  (e) 

gl 

0         for  C  not  in  element  (e) 


-—  for  %   in  element  (e) 

„(e)   _    *e 
g2 


for  %   not  in  element  (e) 
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and  1      is  the  length  of  the  eth  element.   The 


notation 
(i-1) 


in  expression  V.21  means  that  G.  is  the  union  of  g, 
and  g^   .   The  local  (element)  shape  functions  have  the 
following  properties , 


(i)    / 


e   (j)   (m)   _ 


=   0    if   j  ?   m 


0 


(ii)   g<e)(nj)   =   6.. 


if    i  =  j 


if   i  *   j 


Having  defined  the  local  shape  functions,  the  elemental 
matrix  operators  contributing  to  the  global  matrix  operators 
V.16  through  V.20  are, 


-1" 


/  {G.}'  <G.>'dn    % 
o     1   J 


-1 


l  , 

/  (G.}<G.>,dn  |   t 


-1 


-1 


(V.16') 


(V.17») 
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f  {G.}<G.>dn    § 


e 
6 


(V.18») 


-1 


/  {G.  }<G.>'{Y}<G.>dn    4   ^ 


U     2    -><VrVJ 


(V.19') 


/  {G. }dn   fc 


(V.20') 


The  derivations  of  these  operators  are  shown  in  Appendix  E. 

Since  there  is  extensive  coupling  in  the  field  equations 
and  there  are  three  degrees  of  freedom  at  each  nodal  point, 
the  numbering  scheme  for  the  system  matrices  was  important. 
To  minimize  the  bandwidth  of  the  matrices,  the  numbering 
scheme  represented  in  Figure  8  was  used. 

Figure  9  represents  the  matrix  A(t)  of  expression  V.12 
for  any  three  successive  nodal  points.   The  distribution  of 
the  elemental  matrices  for  the  FEM  operators  is  shown,  as 
well  as  the  possible  locations  for  the  coupling  of  the  field 
equations  over  an  element.   In  addition,  Figure  9  shows  the 
bandwidth  that  would  be  observed  for  any  n-1  element  solu- 
tion.  For  this  scheme,  the  bandwidth  is  nine.   If  coupling 
is  not  present,  the  bandwidth  is  seven.   Figure  9  reflects 
the  extensive  coupling  that  is  present  in  the  model.   All  the 
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matrix  elemental  locations  are  filled  except  those  which 
are  shaded  in. 

B.   IMPLEMENTATION  OF  BOUNDARY  CONDITIONS 

Having  formulated  the  system  matrices  for  the  field 
equations ,  treatment  of  the  boundary  conditions  will  now  be 
discussed.   Each  field  equation  is  considered  individually. 
1.   Fiber  Heat  Transfer  Equation 

The  fiber  heat  transfer  boundary  conditions  in  non- 
dimensional  form  from  Appendix  A  are 


n   =   0     X-  -II-  -   h,L(Q-l)  +^(T*-T4)        (V.22) 


"  =  x  Ai  If  "  -h2L(e-1)  "  tt1^-^'     <v-23> 


Since  the  first  term  in  expression  V.13  is 


{Gi)A1<G.>' {0} 


or  in  analogous  form 


Al  3n 


(V.24) 


0 


natural  boundary  conditions  V.22  and  V.23  may  be  directly 
substituted  in  equation  V.24. 
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The  response  dependent  parameters,  h, ,  h2 ,  T  , 
changing  with  time,  are  evaluated  at  the  previous  time 
step.   Thus,  the  boundary  conditions  are  incorporated  in  the 
system  matrices  as  follows, 


(1)  -h,L:   added  to  the  stiffness  matrix  A(t) 

at  location  A,  , 

(2)  h,L  -  ^^(T4  -  T4):   added  to  the  excitation 

1  T    cr    °° 

CO       a 

vector  F(t)  at  location  F, 

(3)  ~n2L:   added  to  the  stiffness  matrix  A(t)  at 

location  A3n_^3n_2 

(4)  h0L  -  ^^(T4  -  T4)  :   added  to  the  excitation 

2  T    g     °° 

CO        ■* 

vector  F(t)  at  location  F-,   0 

2.   Internal  Flow  Heat  Transfer  Equation 

The  non-dimensional  boundary  conditions  presented  in 
Appendix  A  for  the  air  energy  equation  are 


=   0        r   =   1  (V.25) 


n  -  1       M-  -  l~  (V.26) 

<3  n     oil 


The  essential  boundary  condition  at  n  -  0,  r(n  =  0)  =  1  is 

imposed  in  the  Galerkin  equation  as  follows.   The  A9  .  row 

of  the  A(t)  matrix,  the  B0  .  row  and  the  B.  0  column  of  the 
:  2,i  i,2 

B  matrix,  and  the  F?  location  of  the  excitation  vector,  F(t), 
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are  all  set  equal  to  zero.   The  B-  2  location  of  the  B 
matrix  is  then  set  equal  to  one.   The  natural  boundary  con- 
dition at  n  =  1  for  the  air  heat  transfer  equation  was  treated 
in  the  same  manner  as  that  of  the  fiber  heat  transfer  equation, 
Foregoing  the  individual  steps  the  boundary  condition  is 
implemented  by 

h^Lv, 

(1)  -  — ■ :   added  to  the  stiffness  matrix  A(t) 

Al 

at  location  A-,   ,  ~   ~ 
3n-l,3n-2 

h  Lv     aev 

(2)  -=r — =  -  _  :  (T   -  T  )  :   added  to  the  excitation 

A,      I  A,   g     °° 

vector  F(t)  at  location  F^. 

3 .   Oxygen  Transport  Equation 

For  the  oxygen  diffusion  equation,  the  boundary 
conditions  in  non-dimensional  form  from  Appendix  A  are: 


n   -   0  u   it  ■   a,  ($  -  p)  (V.27) 

l  3  n      z 


n  =  1  it  «  o  (V.28) 

3  n 


Since  these  are  both  natural  boundary  conditions,  they  were 
substituted  for  the  first  term  in  expression  V.15  at  n  =  0 
and  n  =  1/  respectively.   Coefficients  w,  and  u)_  are  evalu- 
ated at  the  previous  time  step.   The  boundary  conditions  are 
implemented  by  adding 
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(1)   -<*>-:   to  the  stiffness  matrix  A(t)  at  location 
A3,3 


and 


(2)   w^p:   to  the  excitation  vector  F(t)  at  location 
F3 

This  concludes  the  discussion  for  the  implementation 
of  the  boundary  conditions.   A  word  of  caution  is  in  order. 
After  each  time  step  integration,  the  time-dependent  coeffi- 
cients of  the  boundary  conditions  (i.e. ,  h, ,  h2,    A,  v,  u) 
are  reevaluated.   Before  incorporating  the  updated  coefficients 
into  the  stiffness  matrix  A(t),  the  previous  values  must  be 
subtracted  out.   The  results  of  not  taking  this  into  account 
will  be  obvious . 

The  implicit  system  of  ordinary  differential  equa- 
tions was  integrated  numerically  by  a  modified  implicit-Gear 
method  developed  by  Franke  [2] .   The  time-dependent  coeffi- 
cients of  the  FEM  operators  in  expression  V.13,  V.14,  and 
V.15  were  updated  at  the  previous  time.   The  reaction  rate 
term  (III. 26)  was  also  evaluated  at  the  previous  time,  and 
appears  in  the  excitation  vector  F(t).   In  this  way,  the 
final  system  of  ordinary  differential  equations  was 

A(t)V  +  F(t)   =   B  I 

where  A  and  B  are  (3n  x  3n)  matrices  of  temperature  dependent 
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coefficients.   F  is  a  (3n *  1)  vector  arising  partly  from  the 
reaction  rate  terms  and  partly  from  the  boundary  conditions. 
As  noted,  the  elements  of  F  are  dependent  on  both  temperature 
and  oxygen  concentration. 
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VI.   RESULTS  AND  CONCLUSION 

A.   NUMERICAL  CONSIDERATIONS 

The  input  parameters  which  were  varied  in  the  computer 
analysis  include  (1)  plate  thickness,  L,  (2)  fiber  diameter, 
d,  and  (3)  air  flow  rate,  U^.   Other  parameters  which  could 
be  varied  are  the  ply  thickness,  D,  ambient  temperature,  Ta, 
and  fiber  emissivity,  e.   One  set  of  initial  conditions  was 
used  for  the  analysis.   These  are  shown  in  Figures  10  and  11. 
For  this  initial  effort,  the  actual  initial  conditions  were 
not  of  prime  consideration.   The  selection  of  initial  condi- 
tions of  Figures  10  and  11  will  be  discussed  in  Section  VI. C. 
The  finite  element  program  calculates  the  remaining  system 
parameters  such  as  permeability,  porosity,  pressure  differ- 
ential, and  the  response  variables  as  functions  of  time  and 
position.   Parameters  which  are  functions  of  temperature, 
such  as  k  ,  kq,  k  ,  h.,  p  ,  u  ,  and  y  are  continuously  updated 
during  the  transient  analysis. 

The  preliminary  solution  effort  showed  the  system  of 
equations  to  be  very  stiff  (refer  to  Shampine/Gordon  [24] 
and  Gear  [25]  for  a  discussion  of  "stiffness").   Changing 
the  integration  algorithm  from  a  sixth-order  Runge-Kutta 
(IMSL  subroutine  DVERK)  to  a  modified  implicit-Gear  method 
for  stiff  systems,  developed  by  Franke  [2] ,  resulted  in  a 
significant  reduction  in  CPU  time.   The  computational  effort 
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was  carried  out  on  an  IBM  360/70.   Typical  runs  required  20- 
30  minutes  CPU  time  for  problems  for  which  ignition  occurred, 
and  6-25  minutes  CPU  time  for  extinction  problems.   Ignition 
problems  were  terminated  when  the  temperature  at  any  nodal 
point  exceeded  2200  degrees  Fahrenheit,  or  when  the  graphite 
at  a  nodal  point  was  totally  consumed.   Extinction  problems 
were  carried  out  to  steady  state. 

A  twenty  five  nodal  point  model  (75  o.  d.  e.)  provided 
results  which  differed  less  than  five  percent  from  the  re- 
sults of  a  32  nodal  point  model  (96  o.  d.  e.)  and  was  adopted 
for  all  computer  runs.   For  the  twenty  five  nodal  point 
model,  approximately  275K  bytes  of  core  was  required.   The 
computer  program  which  includes  the  FEM  formulation  and  the 
integration  routine  as  well  as  a  sample  input  file  is  pre- 
sented at  the  end  of  the  appendices. 

The  numerical  solution  produced  satisfactory  results  with 
one  exception  which  will  now  be  discussed.   For  a  particular 

range  of  temperature  and  oxygen  concentration  (approximately 

3 
1500  degrees  Fahrenheit  and  .001  lbm/ft  ,  respectively), 

there  is  a  significant  increase  in  reaction  rate.   This 

accelerated  reaction  rate  produced  a  negative/positive 

oscillation  for  the  nodal  values  of  oxygen  concentration. 

The  oscillation  occurred  in  a  region  of  the  plate  where  the 

oxygen  appeared  to  be  totally  consumed.   As  discussed  by 

Frank-Kamenetskii  [3],  nth  order  reactions  may  yield  zero 

values  for  the  oxygen  concentration  and  gradient  within  the 
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medium.   The  cause  of  the  oscillation  was  as  follows. 
Although  the  reaction  rate  is  a  function  of  temperature  and 
oxygen  concentration,  numerically,  it  is  treated  as  constant 
during  an  integration  interval.   As  a  result,  the  concentra- 
tion in  the  region  of  high  reaction  rate  becomes  negative. 
The  reaction  rate  was  updated  for  the  next  time  interval  by 
setting  negative  concentrations  to  zero.   During  this  time 
interval,  for  the  region  of  zero  concentration,  there  is  no 
reaction.   Without  consumption,  the  oxygen  concentration  will 
increase  and  the  cycle  repeats  itself.   The  oscillation  was 
aggravated  by  large  plate  thickness  and  low  permeability. 
The  instability  occurring  for  the  high  reaction  rate  may  be 
corrected  by  (1)  numerically  integrating  the  rate  term  in 
the  time  domain,  (2)  iterating  until  convergence  is  obtained 
between  consecutive  values  of  concentration,  or  (3)  decreasing 
the  time  step  and  the  length  of  the  elements.   As  a  matter 
of  expediency,  measure  (3)  was  used  to  minimize  the  insta- 
bility.  The  fiber  and  air  temperatures  were  essentially 
unaffected  by  the  oscillations  in  the  oxygen  concentration. 

B.   RESULTS  AND  OBSERVATIONS 

Table  1  presents  the  results  of  fourteen  problems.   In 
Table  1,  the  parameters  which  depend  upon  temperature  (u , 
Re.,  and  h.)  are  evaluated  at  ambient  conditions  for  com- 
parison  purposes  only.   These  parameters,  in  fact,  will  vary 
during  the  course  of  the  transient  analysis.   In  all  cases, 
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the  initial  condition  on  graphite  temperature  was  1050 
degrees  Fahrenheit. 

Figures  10,  11,  and  12  show  the  transient  behavior  of 
Case  1.   For  the  initial  conditions  shown,  ignition  occurred 
in  approximately  eight  seconds.   The  reaction  rate  and  oxygen 
concentration  responses  for  Case  1  are  typical  of  the  problems 
for  which  ignition  resulted.   Figure  13  shows  that  the  air 
and  fiber  temperatures  for  Case  1  do  not  differ  significantly. 
This  behavior  is  typical  for  most  problems.   However,  for 
plates  with  high  porosity  (greater  than  .8)  and  subjected  to 
relatively  high  flow  rates,  significant  differences  in  the 
air  and  fiber  temperatures  do  occur  (see  Figure  14  for  Case 
9) .   In  a  previous  effort  by  Vatikiotis  and  Salinas  [26] , 
linearly  varying  initial  conditions  were  investigated.   Although 
the  reaction  rate  used  in  that  investigation  differed  from 
the  one  used  in  the  present  investigation,  the  overall  re- 
sults remain  valid.   The  main  observation  was  that  ignition 
is  less  likely  to  occur  for  thin  plates  with  the  ambient  air 
entering  the  hotter  plate  surface.   Since  temperature  con- 
trols the  reaction  in  the  kinetic  regime,  cooler  air  entering 
the  hotter  plate  surface  enhances  the  heat  loss.   This  reduces 
the  fiber  temperature,  thereby  decreasing  the  reaction  rate. 
In  effect,  the  cooler  air  enters  the  plate  where  it  is  most 
needed. 

In  the  present  analysis,  observations  were  made  by  vary- 
ine  one  input  parameter  (i.e. ,  L  or  d  or  U  )  while  keeping 
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all  others  fixed.   The  effects  of  each  parameter  on  the 
behavior  will  be  discussed  individually. 
1.   Effects  of  Exterior  Velocity 

As  shown  in  Figures  15,  16  and  17  (cases  4,  1  and 
2,  respectively)  for  an  initial  condition  of  1050  degrees 
Fahrenheit,  three  distinct  regions  of  combustion  were  ob- 
served (extinction-ignition-extinction) .   These  became  appar- 
ent as  U  was  increased  from  low  velocities  (10  knots)  to 

00 

higher  velocities  (120  knots) .   Vulis  [4]  discusses  an  experi- 
ment of  a  heated  carbon  rod  with  an  air  jet  impinging  upon 
it.   For  a  certain  range  of  flow  velocities,  ignition  was 
observed.   However,  extinction  did  occur  for  velocities  less 
than  and  greater  than  the  velocity  for  which  ignition  occurred. 
The  behavior  of  the  carbon  rod  and  that  of  the  graphite  mat 
is  similar.   This  behavior  may  be  explained  by  the  Semenov 
model  of  Figure  1.   For  high  Uw,  the  internal  heat  transfer 
coefficient,  h.,  will  be  large  as  a  result  of  an  increase 
in  the  pore  velocity.   Since  the  slopes  of  the  heat  loss 
lines,  q  ,  increase  with  increasing  h.  for  a  constant  T  , 
the  graphite  temperature  decreases  from  the  critical  ignition 
point  I.   This  causes  extinction  for  the  higher  exterior  flow 
velocities,  U^.   In  contrast,  for  low  U  ,  the  slope  of  q  will 
be  small  and  ignition  is  more  likely  to  occur.   However,  the 

heat  generation  curve,  q  ,  changes  position  due  to  a  decrease 

y 

in  oxygen  concentration.   This  reduction  in  oxygen  concentra- 
tion is  caused  by  the  lower  induced  flow  rate  through  the 
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plate.   The  overall  shape  of  q   is  flattened  such  that  the 
q„  line  crosses  q  behind  the  new  critical  point  I.   Thus 
extinction  is  observed  for  this  region  of  low  U  . 
2.   Effects  of  Fiber  Diameter 

Varying  the  fiber  diameter  with  the  remaining  input 
parameters  fixed,  also  yields  three  regions  of  combustion 
(extinction-ignition-extinction)  for  a  fixed  initial  tempera- 
ture.  The  results  are  shown  in  Figures  16,  18  and  19  (Cases 
1,  9  and  10,  respectively).   The  fiber  diameter  and  ply 
thickness  determine  the  permeability  of  the  plate.   Similar 
to  U^,  permeability  affects  the  pore  velocity  (i.e. ,  low 
permeability  -»■  low  velocity)  .   In  turn,  this  affects  the 
convective  heat  transfer  and  the  amount  of  oxygen  entering 
the  plate.   These  three  regions  of  combustion  are  explained 
by  the  Semenov  model  for  the  same  reasons  as  those  of  varying 
U  . 

00 

Porosity  is  a  measure  of  void  space  and  denotes  the 
space  available  for  oxygen.   Porosity  is  associated  with  the 
internal  geometry  of  the  medium,  as  is  permeability.   However, 
they  are  basically  different  parameters.   Permeability 
(hydraulic  conductivity)  is  associated  with  convective  air 
flow  through  the  medium.   Thus,  an  n-fold  increase  in  the 

fiber  diameter  and  ply  thickness  will  not  affect  a  change 

2 
in  porosity,  but  will  affect  an  n  order  change  in  the  permea- 
bility.  For  an  n  greater  than  one,  the  result  is  a  decrease 
in  pore  velocity  without  a  change  in  porosity.   Porosity 
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determines  the  maximum  oxygen  concentration  per  unit  volume 
of  fibrous  mat.   This  indirectly  affects  the  reaction  rate. 
Thus,  for  a  given  u  ,  the  lower  the  porosity,  the  less  oxygen 
there  is  for  combustion.   It  appears  that  high  porosity 
would  enhance  combustion  by  the  presence  of  more  oxygen. 
However,  higher  porosities  also  provide  more  fluid  per  unit 
volume  resulting  in  greater  heat  transfer.   This  is  observed 
in  Case  9,  Figure  18. 

3.   Effects  of  Plate  Thickness,  L 

For  the  given  initial  temperature  of  1050  degrees 
Fahrenheit,  varying  the  plate  thickness  yielded  two  regions 
of  combustion  (extinction-ignition) .   Extinction  was  observed 
for  thin  plates  (Case  2,  Figure  17) .   As  shown  by  Figures  20 
and  21  (Cases  7  and  11,  respectively),  ignition  was  observed 
when  plate  thickness  was  increased.   Plate  thickness  also 
affects  the  pore  velocity.   Decreasing  L,  increases  the 
pressure  gradient  across  the  plate,  thus  increasing  the  pore 
velocity.   The  result  is  that  extinction  is  more  likely  for 
thin  plates  because  of  the  enhanced  convection  heat  trans- 
fer resulting  from  the  increase  in  pore  velocity.   Whereas 
pore  velocity  is  inversely  proportional  to  L,  it  is  an  exponen- 
tial function  of  U^.   For  the  range  of  parameters  in  this 
investigation  (i.e. ,  .1"  <  L  <  3.",  and  10  knots  <  U^  <  100 
knots) ,  an  n-fold  change  in  Uot  will  produce  a  greater  change 
in  u  ,  then  will  an  n-fold  change  in  L.   Compare  Cases  1  and 
2  to  Cases  7  and  8. 
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A  characteristic  feature  of  ignition  for  thick 
plates  was  observed.   As  the  thickness  increased,  the  igni- 
tion region  (the  spatial  location)  moved  closer  to  the 
entrance  surface  of  the  ambient  air.   Compare  Figures  20 
and  21. 

C.   DISCUSSION 

At  this  point,  discussion  is  made  of  ignition  tempera- 
tures for  porous  graphite  plates.   The  ignition  tempera- 
ture for  Case  1  was  1050  degrees  Fahreheit.   This  tempera- 
ture was  obtained  by  varying  the  initial  temperature  until 
ignition  occurred.   In  all  subsequent  cases,  1050  degrees 
Fahrenheit  was  adopted  as  the  initial  condition,  and  we 
observed  whether  extinction  or  ignition  occurred.   The  igni- 
tion temperatures  given  here  are  the  result  of  taking  n  =  -s- 
in  the  nth  order  reaction  rate.   Significantly  different 
ignition  temperatures  will  result  for  other  n.   For  example, 
for  n  =  1  (1st  order  reaction) ,  the  ignition  temperature  for 
Case  1  was  approximately  1600  degrees  Fahrenheit.   As  Frank- 
Kamenetskii  [3]  has  suggested  values  of  n  between  1/3  and 
2/3,  the  average  value  of  1/2  was  adopted  for  this  analysis. 
The  ignition  temperatures  obtained  using  n  =  -j   agree  favorably 
with  the  experiments  of  Fontenot  [1] .   For  fiber  combustion, 
ignition  temperatures  for  porous  graphite  plates  are  on  the 
order  of  1100  degrees  Fahrenheit. 

Consider  the  results  of  varying  U^.   Extinction  occurred 
for  Case  2  and  Case  4;  thus,  the  ignition  temperature 
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for  each  case  is  greater  than  1050  degrees  Fahrenheit. 
Since  ignition  does  occur  at  1050  degrees  Fahrenheit  for 
intermediate  values  of  U^  (Case  1,  Case  7),  ignition  tem- 
perature is  not  a  monotonic  function  of  U  .   The  precise 
determination  of  ignition  temperature  for  a  specific  value 
of  air  velocity  can  be  obtained  by  a  large  number  of  com- 
puter runs.   A  similar  discussion  can  be  made  for  the  igni- 
tion temperatures  dependence  on  fiber  diameter.   The  general 
shape  of  these  functions  is  represented  by  the  shaded  region 
on  Figure  22.   These  suppositions  may  only  be  valid  in  the 
limited  range  of  parameters  associated  with  the  present 
investigation . 

It  is  interesting  to  observe  that  changing  the  initial 
condition  of  the  oxygen  concentration  only  affects  the  early 
part  of  the  problem.   The  response  of  the  concentration  is 
fast  compared  to  that  of  temperature.   After  approximately 
one  second,  the  transient  behavior  is  the  same  regardless 
of  the  initial  condition  on  the  oxygen  concentration. 

A  brief  description  of  the  Semenov  model  was  given  in 
Section  II. A.   That  discussion  was  for  a  lumped  parameter 
model.   The  actual  combustion  process  is  more  complex  since 
there  are  spatial  variations  in  temperature  and  in  oxygen 
concentration.   The  results  of  this  analysis  show  that  both 
kinetic  and  diffusion  regimes  may  exist  simultaneously  in  a 
porous  medium.   Hence,  it  is  not  reasonable  to  restrict  an 
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analysis  of  combustion  in  a  porous  medium  to  the  kinetic 
regime. 

D.   CONCLUDING  REMARKS 

The  model  has  provided  considerable  insight  as  to  the 
behavior  of  fibrous  composites  subject  to  combustion. 
Improvements  may  be  realized  by  extending  the  model  to 
include: 

(1)  The  combustion  of  the  epoxy  matrix  (i.e./  start 
problem  with  first  stage  of  combustion  process) . 

(2)  The  change  in  geometry  of  the  porous  medium  due 
to  combustion  (affects  p,  m  and  z) . 

(3)  A  more  complex  reaction  (i.e. ,  secondary  combustion 
of  gaseous  by-products,  such  as  carbon  monoxide). 

(4)  The  effect  of  gaseous  by-products  on  fluid  properties 

(5)  Generalization  to  two  and  three  dimensional  models 
(results  in  anisotropic  properties) . 

The  results  show  that  certain  properties,  such  as  plate 
thickness,  permeability  and  porosity,  have  significant  effects 
on  combustion.   Therefore,  one  could  select  these  design 
parameters  to  improve  the  flame  resistance  of  composite 
materials.   This  would  be  beneficial  to  the  survivability 
of  a  composite  structure.   A  comprehensive  set  of  analyses 
could  be  undertaken  to  couple  the  combustion  problem  to  the 
strength  problem.   In  this  way,  one  could  achieve  a  design 
which  maximizes  strength  and  minimizes  combustion. 
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Another  use  of  the  present  model  would  be  in  the  study 
of  energy  systems  utilizing  particulate  fuels  (i.e. ,  coal, 
biomass) .   In  this  case,  the  design  parameters  could  be 
selected  to  maximize  the  performance  of  the  system. 
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APPENDIX  A 
NONDIMENSIONALIZATION  OF  FIELD  EQUATIONS 

The  field  equations  and  boundary  conditions  are: 

3T  h.  Z  3T 

f-[(k  +k    )-r-S]  -  rr~r(T  -T, )  +  AHr(T    ,»)    =    P    c-^        (III. 28) 
3x        g      r     3x  (l-p)      g      a  g  g  g  3t 

3T  3T  h.z  3T 

k<kaTT'  "  ™caTT  +  ^-(VTa»    =   pacalF  (III"37> 

k<4£>    -  upH  -  Sjhk>*   -    <|)r(Tg,*)      -     ||  (III. 46) 


at   x  =   0 


3  T 

(k+k)-J      =      h1  (T    -T    )    +    ae(T*-T*)  (IV. 1) 

gr3x  10=°  goo 


T        =      T  (IV. 3) 

a  °° 


B   H     =      Up(*-P^)  (IV-5) 


at   x   =  L 


3  T 

(k   +k    )-r-2      =      -h    (T   -T    )     -    oe(T*-T*)  (IV. 2) 

grdx  zg°°  g°° 


3T  3T 

__*      =      _£ 

3x  3x 


(IV. 4) 
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|4   =   0  (IV. 6) 

j  x 


The  nondimensional  variables  are  defined  as: 

0   =   T  /T     nondimensional  fiber  temperature 
y 

r   =   T  /T     nondimensional  air  temperature 
a  °° 

*   =   ^/^oo     nondimensional  oxygen  concentration 

n   =  x/L      nondimensional  distance 

The  time  variable,  t,  will  not  be  nondimensionalized. 

Using  the  temperature  of  the  fiber,  T  ,  to  demonstrate 
the  technique  of  transformation,  we  have, 


T    =   T  0 
g       oo 


It  follows  that, 


3T 

— 3.      =   T   1© 

3x  <=°  3X 


92T  „2„ 

1  =  T  14 

2  <*>    2 

3x  3x 


Using  the  chain  rule  to  transform  x  to  n ,  gives 


3T 

g  =  T  i©  in 

3X        oo  3n  3x 
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Since  n  =  x/L,  the  partial  derivatives  become 


and 


Therefore, 


3x      L  3n 


32T       T  _2_  , 

g  _    «>  3 0  _3_Q. 

^2      ~L*  I    2  3x 

3x  3n 


92T  T   .2. 

g      _°°_  3 0 

2  ~    2    2 

3xz  L   3n 


Similar  transformations  for  T   and  4>  yield, 

a 

a      _        <»    3  r  3$      _      Toa    3  $ 

~3x~  T"  Tt\  3~x      ~      ~L~   3n 

32T       T   a_  .2.      <J>   .2. 

„  2       .2  9ri  „  2      _2„2 

3x       L  3x      L   3n 

Substituting  these  relations  into  the  field  equations  and 
removing  the  constants  from  the  differential  operators  gives, 

T  h .  z 

z   3n   g  r  3n    (l-p)  g        g  got 


T   a     or    mc  T       h.  z 

_2  3n   a3n      L    3n    p  a  a3t 

Li 
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^2  37(B— }  "   L   37T  -  -L-(3T^3TT*  "  (f)r(Tg'*}   "   ♦-at 


Upon  rearrangement,  the  equations  become, 

2 

3      r/i,    j.i      ^Qi  hiZL     ta     n^  AHL2     ,_       aN  t23G 

^WalT1  "  TI^?y(0"r)  +  "tT^V40     =    pgcgL  at 
3    ,,    3i\      *     T3r       hizL   #ri   _.  T23r 

t—  (k    -r— )    -  mc    Lt— -   +    (9-r)        =       P     C    L    t-z- 

3na3n  a    3n  p  a   a      3t 

3     ,n3*.  T3$       _    T,3Up,3rA  /1>L       ,m       ,»  T2    3$ 

ail        d  n  p     on  °°       3T       3  r|  r     <p  g  3t 


a 


Letting, 


h   zL2  . .,. : 


X.      =      (k     +k      )                   X~     =     -ynr r-                 A  _     =     — = A  „      =     p      C     L 

1             g      r                2         (1-p)              3           T^  4         Mg   g 

h.  zL  ~ 

i                                         t                             i  t  2 

v,       =      k                     v_    =  mc   L              v_    =   — - v  .    =    p    c   L 

1               a                   2a                 3P  4          a  a 

3u  ,       2 

u>,    =   B                          con    =    u  L                u-    =    T   L7r=2.  u      =   -  — 

1                                       2p                      3           °°    3T  4         f    <j> 

a  °° 


lllr       =      L 


2 


the   field  equations   may  be  written   as, 


k(Xl  $    "    X2(9"r)    +    V(Tg'*>      =      X4   H  (A-1' 


1^1   &    —  2|+    «3(9-r»       "      -4   li  (A-2> 
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37<u,i  aT1   -  u2  an  "  "3   3?  -  V'V*'     =     "5  at     <A'3) 


Noting  that  the  coefficients  \,    v,  co  are  response  dependent 
variables ,  equations  A.l,  A. 2,  and  A. 3  are  considered  in 
their  final  form  in  Section  III. 

Substituting  the  nondimensional  variables  into  the 
boundary  conditions  gives, 


and 


r  =  l  x  =  o 


wl^  =   <V$-p)  X  =  ° 


x  |£  =  .hl(e-i,  -  ¥&!**  -  ii) 

m       3 


9T   =   36 

3n     3n 


x  =  L 


r—   =   0  X  =  L 

3n 


The  dimensional  form  of  the  absolute  temperature  will  be 
retained  for  convenience. 
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APPENDIX  B 
POLYNOMIAL  APPROXIMATIONS  OF  THERMAL  PROPERTIES 

Relations  giving  the  viscosity  and  thermal  conductivities 
of  air  for  varying  temperature  were  required.   A  simple  way 
to  obtain  values  for  these  properties  is  to  fit  empirical 
data  with  2nd  order  Lagrange  polynomials. 

The  general  form  of  the  2nd  order  Lagrange  polynomials 
for  thermal  conductivity  and  viscosity  are 


(T  -T   )  (T  -TO  (T  -T  -.)  (T  -T  ,) 

k        a   a    a   a3    ,        a   a3   a   al 


(Tal-Ta2}  <Tal-Ta3>   al    (Ta2-Ta3) (Ta2-Tal}   *2 


(T  -T  ,) (T  -T  9) 
+    a        a   a      V  (n    i  "\ 

(T  ,-T  ,)  (T   -T  .)   a3  [i5ml) 

a>j   a_L   a  j   a^ 


(VTa2}  (VTa3}         (VTa3}  (VTal} 
"  "   (Tal-Ta2} (Tal"Ta3)  ^         (Ta2-Ta3>  «Ta2"Tal)  "2 


(T  -T  .) (T  -T  0) 
a   al   a   a2 

(T  _-T  ,)  (T  ,-T  75  U3 
a3   al   a3   a2 


(B.2) 


Choosing  three  points  from  a  range  of  temperatures  that 

would  be  representative  of  those  observed  during  the  analysis, 

the  corresponding  values  of  the  properties  are: 
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(subscript) 

Temp 

(F) 

1 

0 

2 

700 

3 

1500 

ka 

,    Btu    >  ,  lbm  s 

llbm-hr-ft;  lft-hrJ 


.0131  .0394 

.0284  .0765 

.0423  .107 


Substituting  these  values  into  expressions  B.l  and  B.2 
gives, 


(T  -700) (T  -1500)  (T  -1500) (T  -0) 

ka   =    (0-700)  (0-1500)  (-0131)  +  (700-1500)  (700-0)  ('0284) 


(T  -0) (T  -700) 

+  5 £ (    0423} 

(1500-0)  (1500-700)  UU^J; 

(Ta-700) (T  -1500)  (T  -1500) (T  -0) 

y   =   (0-700)  (0-1500)   (-0394)  +  (1500-0)  (1500-700)  ('0765) 


(T  -0) (T  -700) 

(.107) 


(1500-0) (1500-700) 
These  expressions  reduce  to. 


k   =   -3.232 x   10"9T2  +  2.412  xl0"5T   +  .0131       J  J 


'a  '  ~  a  "  "  a  lbm-hr-ft 


y      =      -1.041x10    8T2   +    6.029x10    5T      +    .0394 


a  '    "'  a  ft-hr 


82 


Both  expressions  give  property  values  which  are  within 
five  percent  of  the  data  for  temperatures  to  2000  degrees 
Fahrenheit. 
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APPENDIX  C 
TRANSFORMATION  OF  THE  REACTION  RATE  TEEM 

The  expression  for  the  oxidation  rate  of  graphite  fibers 
taken  from  Parker  and  Hottel  [15]  is 


g 
9.55x10   _       ,-44000,        gm         ,„    ,. 

rrr   =   TTTo pno  exP( : )  T (C.l) 

9       T1/2     °2     R  T  .  cm2-sec 
gk                gk 


where  R  is  the  universal  gas  constant  (1.958  cal/gmole-K) 
and  PQ2  is  the  partial  pressure  of  oxygen  (atm.).   The 
partial  pressure  may  be  represented  in  terms  of  the  tem- 
perature, T,  and  oxygen  concentration,  <j> ,  by  the  Ideal  Gas 
law  in  the  form  given  by  Kanury  [16] 


M 
P02   =   <M^>  RA  T  *  <C-2) 


where  M  /M02  =  • 9  is  the  ratio  of  the  molecular  weights  of 
air  and  oxygen;  R   (=  53.34  ft.-lbf/lbm  R)  is  the  gas 
constant  for  air.   Substituting  M  /Mn2  into  equation  C.2 
and  converting  P_2  to  atmospheric  units  the  expression 
becomes 


P_n   =   4.25  x  10~4  R-  T  <j>     atm. 
02  A    r 
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2  2 

Converting  r   to  lbm/ft  -hr  from  gm/cm  -sec,  equation  C.l 

becomes 


7.04  x  1Q1Q         ,-44000,        lbm 

r~  T~~ 70        ^nO   eXP  I   D   T<     ' 


9  *5r  R  V        ft2-hr 


Substituting  in  for  Pn~  and  converting  the  absolute  tem- 
peratures from  Kelvin  to  Rankine,  the  rate  expression 
becomes, 


,  ft,,   in7  RA  T  *     ,-39883,      lbm 
r   =   4.014  x  10   y-p?   exp  ( — - ) 


9      "*'     '    T1/2       T       ft2-hr 

g       g 


The  absolute  temperature  at  the  fiber  surfaces  is  T  .   Thus, 
T  in  the  numerator  becomes  T  ,  and  the  expression  becomes 

g 


r    =   4.014  x  107  R   T1/2  <j>  exp(  39883) 

g  A  g    y    *        T 

g 

3 
To  transform  the  reaction  rate  to  lbm/ft  -hr,  y,  the  ratio 

of  a  fibers  surface  area  to  its  volume  is  placed  in  the 

numerator  of  r  .  The  expression  for  y  is, 


1  TT  d       2  ,_,  -,, 

y  =  it?  -  d  (c-3) 

4* 

where  d  is  the  fiber  diameter,  and  the  y  factor  accounts  for 
the  uncertainty  of  the  actual  amount  of  fiber  surface  after 
the  epoxy  has  burned  away.   The  reaction  rate  of  the  graphite 
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fibers  per  unit  volume  is  expressed  as 


a      ni a         in7     n    m^/2  ,"39883,        lbltl 

r   =   4.014  x  10   y  R   T  '  $    exp  ( — )  — = 

A      g  T  ft  -hr 

g 


To  obtain  the  heat  liberated  per  lbm  of  fiber,  r  is 
multiplied  by  AH,  the  enthalpy  of  formation  of  carbon  dioxide 
per  unit  mass  of  graphite  consumed.   This  value  is  approxi- 
mately 14,085  BTU  per  pound-mass  of  fiber  and  was  obtained 
from  Kanury  [16] 


Heat  generation  rate   =   AH  r(T  ,  4>) 


To  obtain  the  oxygen  consumption  rate,  r  is  multiplied 
by  the  inverse  of  the  fuel-air  stoichiometric  ratio, 


Oxygen  consumption  rate   =  (-?)    r(T  ,<£). 
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APPENDIX  D 
JUSTIFICATION  OF  THE  DANKWERTS f  BOUNDARY  CONDITIONS 

The  Dankwerts '  boundary  conditions  for  the  oxygen  diffu- 
sion equation  are: 


dx 


u 
B  v 


"  P4>  > 


x  =  0 


d» 
dx 


=   0 


x  =  L 


Bischoff  [23]  presents  a  discussion  of  these  equations.   A 
brief  summary  is  provided  for  completeness. 

Figure  10  shows  the  one-dimensional  oxygen  transport 
region  considered  in  this  analysis. 


x  =  0 


x  =  L 


Flow 


Entrance 
Section 


4>1/    B^    u1 


Porous 
Plate 


,  B,  u 


Exit 
Section 


♦2,  B2,  u2 


The  concentration  of  oxygen,  <j>,  for  each  region  is  distin- 
guished by  subscripts;  B  and  u  are  the  diffusivity  and  the 
velocity,  respectively,  for  each  region.   The  oxygen  diffu- 
sion equations  for  each  section  are  as  follows, 


87 


Bl  d2»l    d* 

ul    dx2         dx 


=      0 


x    <    0 


(D.l) 


_B_  d_*    _    d£   _        (  }       u       0 

Up    dl7         dX 


0    <    x    <    L  (D.2) 


^d2*        ^2 
u2   dx1"         dx 


=      0 


x    >    L 


(D.3) 


with   the   general   boundary   conditions, 


(-oc)  =  A 


(D.4) 


p<J>1(0    )       =      <j>((T) 


(D.5) 


p^(o  )   -  p  —      dx 

1 


Bld*l<°>     =     ,<„♦,   -Adiio!! 


u  dx 

P 


(D.6) 


p<j>(L    )       =      4>2(L    ) 


(D.7) 


P*(L  )-p- 33^—    -     *2(L  )    -  —  — 3J 

P  2 


(D.8) 


<j)2(+oo)      =      finite 


(D.9) 
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For  convenience,  the  parameters  u.,  B.  will  be  treated  as 
constant.   In  addition,  the  porosity,  p,  has  been  included 
to  account  for  the  discontinuity  that  exist  for  the  concen- 
tration.  The  discontinuity  is  caused  by  the  oxygen  from  a 
totally  gaseous  environment  to  one  that  is  partly  occupied 
by  solid. 

An  analytical  closed-form  solution  for  this  set  of 
equations  is  not  possible  because  of  the  nonlinearity  of 
equation  D.2.   However,  the  solutions  of  equations  D.l  and 
D.3  are, 

<J>,   =   A,  +  A2  expCUj^x/B,)  (D.10) 

*2   =   A3  +  A4  exP(u2x/B2)  (D.ll) 

Applying  boundary  conditions  D.4  and  D.9,  the  following  re- 
sults are  obtained, 


Al 


A4   =   0 


The   solutions   D.10   and  D.ll   become 


<j>,      =      <j>      +  A0   exp(u.x/Bn)  (D.12) 

1  °°  z  11 
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<J>2   =   A3  (D.13) 


It  would  be  necessary  to  have  the  solution  for  the  non- 
linear equation  D.2  to  solve  for  A-  and  A-.   However,  the 
constants  need  not  be  known  to  continue  with  the  analysis. 
From  equation  D.12, 


4>1(0    )       =      *m   +   A2 


and, 


d<j>,  (0    )  u 

i =  x    a 

dx  B1      2 


Substituting   these    into   equation   D.9   gives 


p«      +  pA,  -  p  !i(!ji  A,)      =      *(0+)-A^i^I 
xry°°        r   2  u,    B,      2  Y  u  dx 

11  P 


Cancelling  terms  yields, 


wn+»    B  d^>(0  ) 

^T°°         TV/      u       ^x 

P 


Rearranging,  the  Dankwerts '  boundary  condition  at  x  =  0  is 


«±  .   J£U  -  P*J  x  =  0 
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From  equation  D.13,  and  noting  that  A-  is  a  constant,  we 
have 


d4>9 

— £  =   0 
dx 


and 


d<j>9(L+) 

— 3 "   ° 

dx 


Substituting  these  expressions  and  equation  D.7  into  equation 
D.8, 


POU")  -  p  A  ^-1   =   p*(lT)  -  ^(0) 
P  2 


Upon  cancelling  terms,  the  second  Dankwets *  boundary  condi- 
tion becomes, 


Si     =  0  x  =  L 

dx 


An  important  underlying  consideration  for  using  the 
Dankwerts '  boundary  conditions  is  that  it  simplifies  the 
analysis  since  equation  D.2  may  be  solved  independently 
without  haveing  to  consider  entrance  and  exit  regions. 
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APPENDIX  E 
DERIVATION  OF  THE  FEM  OPERATORS 

In  the  section  on  the  finite  element  formulation,  the 
following  five  differential  operators  were  identified, 


1 
/    {G!}<G'.>dn  (E.l) 

nil  * 


1 

/    {G.}<G'.>dn  (E.2) 

n      1        3 


1 
/    {G.}<G.>dn  (E.3) 

«  1  1 


1 

/  {G. }<G!>{y}<G.>dn  (E.4) 

n   l    3       3 


1 
/  (Gi}dn  (E.5) 


where  the  G.  are  the  global  basis  functions.   These  opera- 
tors are  constructed  on  the  element  level  by  introducing  the 
corresponding  element  basis  functions,  g. .   The  global  and 
element  basis  functions  are  related  by, 
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r        _      „(i-D  &     „(i) 
Gi  -  g1  9  g2 


where  g.  and  g~  are  defined  by 


gl 


g2 


(1  -  -~)    for  ?  in  element  (e; 
*e 


for  c  not  in  element  (e) 


-p-         for  5  in  element  (e) 


for  %   not  in  element  (e) 


and  l      is  the  length  of  the  (e)th  element. 

The  derivation  of  the  local  elemental  matrices  according 
to  the  Galerkin  method  for  the  global  operations  E.l  through 
E.5  proceeds  as  follows: 
For  operator  E.l, 

global  local 


i  1b  |»i. 

/    {G-  }<G!>dn  %        U[    /        '  <g-   g'>dU 


Noting  that, 
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the  elemental  matrix  becomes 


<f>2    -<f>2 

e  e 


-(f)2  <f>2 

e  e 


dC      I 


-1  1 


-1       1 


For  operator  E.2, 


global 


local 


/  {Gi}<G!>dn    |     U[  / 


e   ^1 


<g^  g£>dc] 


Substituting  in  the  local  shape  functions  gives , 


(i  -f-) 

e 


<f> 

e 


e   e 


Carrying  out  the  operations  gives 
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(-  -i-  +  -L) 

e    e 


-(-**> 


e    e 


{*-*) 


21' 


21 


«    I     I 


r  -l    1 1 


-l 


For  operator  E.3, 


global 


local 


/  {Gi)<G.>dn    }    U[  / 


£e  I*!] 


<gx  g2>d£] 


Substituting  in  for  the  local  shape  functions 


(i  -f-) 

e 


e 


<(1   - 


e        e 


the   elemental  matrix  becomes 


r  (i-  2 


-£-  +  *-) 

£2 

e      jT 

e 


*  „2; 


d-2  X+  i^j 

*e        jT 
e 


dc     I 


r  i 
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For  operator  E.4, 


global 


local 


J-  . 

/  {G. }<G!>  {Y}<G.>dn    W  U{  / 


£   g- 


e   *1 

■i»2 


«31<32>d^ 


y 


Multiplying  out  the  matrices,  the  expression  becomes, 


r(gigigi*i-i +  gigig2 V  (gigiVi-i +  gig2g2 V 


L(glglg2  Vl  +  glg2g2V        (gig2g2Vl  +  g2g2g2V  J 


d? 


Integrating  each  term  separately, 


/      g^g-ids    = 

0 


£72 

e        £ 

e 


e 


1 

3 


/    ^g^^dc    =      J 


e  e 


1 
3 


fQ  gigig2d? 


e  2  , 

0  %l         i*  6 

e  e 
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/    sig2g2d? 

0 


e 
0 


4 


e 


1 
6 


L   glglg2d^   =   -  6 


'0  glg2g2dC   "   I 


JQ  gig2g2dc  = 


0     <T 


/  g-.g;>g?<H  = 
o 


^3  cU 


0   £ 


Substituting  the  values  into  the  expression 


<"  I  Vi  +  7  V 


("  I  Vi  +  I  V 


("  I  Vl  +  k   V 


("  I  Ti-1  +  7  V- 


factoring  out  a  -  1/3,  the  local  matrix  for  the  global 
operator  becomes 
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1  -, 

/    {G. }<G.>{y}<G.>dn  t  "   T 


(Vi-v    (^i — i)' 


'    2    -'    "i-l'V- 


For  operator  E.5, 


global 


local 


1  e  ^1 

/  {Gi)dn    |     u[  J    . 


Substituting  in  the  local  shape  function,  and  integrating, 
the  expression  becomes 


e 


l 


«     I 


This  last  operator  is  used  for  the  excitation  vector  as 
described  in  the  FEM  formulation. 
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FIGURE  2  Idealized  geometry  of  a  fibrous  graphite 
plate. 
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into  respective  volumes  of  fiber  and  air. 
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FIGURE  7  Linear  shape  functions  for  the  Galerkin  formulation 
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FIGURE  9  Shows  coupling  in  system  of  p.  d.  e. 
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FIGURE  22  Shaded  region  represents  shape  of  ignition 
temperature  curve  as  a  function  of  U^or  d. 
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